Tema 6 Inferencia estadística

INFERENCIA ESTADÍSTICA

La inferencia estadística es el conjunto de técnicas estadísticas que nos permiten llegar a conclusiones sobre una población a partir de los datos de una muestra de dicha población. La estadística por tanto parte de los datos para extraer unas conclusiones, es un método deductivo. Para dar soporte teórico a la estadística existe el cálculo de probabilidades, que lo que proporciona es un método inductivo, es decir desde una regla es capaz de llegar a la misma conclusión. Tanto para la estadística como para la probabilidad, existen dos corrientes:

Frecuentista: Se basa en la repetición continuada de los hechos. Cuanto más se repita un experimento, la probabilidad de que ocurra un suceso será más regular y su frecuencia relativa se acercará mucho a la probabilidad del suceso.

USAREMOS ESTE MÉTODO CUANDO NO SEPAMOS ANTES LAS PROBABILIDADES DE CADA SUCESO, POR EJEMPLO QUE UNA MONEDA ESTÉ TRUCADA Y NO SEPAMOS CÓMO, SÓLO POREMOS SABER SU PARCIALIDAD REPITIENDO, POR ELLO SE LLAMA CONOCIMIENTO A POSTERIORI

La interpretación clásica, mayoritaria por lo menos hasta ahora, define la probabilidad en términos de experimentación. Si repites un experimento un número infinito de veces y compruebas que en 350 de cada 1.000 ocasiones se ha producido un determinado resultado, un frecuentista diría que la probabilidad de ese resultado es del 35%. Basándose en esta definición, un frecuentista afirma que es posible asociar a cada evento una probabilidad de obtener un valor VERDADERO del mismo.

La aproximación clásica se basa por lo tanto en estudiar la probabilidad real de las cosas, tratando de determinar hasta qué punto una medición realizada sobre un conjunto de experimentos se aproxima a la probabilidad real que subyace. Es por ello que un frecuentista definiría probabilidad como una expresión matemática que predice en qué medida es verosímil que ciertos eventos ocurran basándose en el patrón observado hasta este momento.

El enfoque frecuentista considera que los datos son aleatorios debido a la naturaleza del proceso que los genera y que cada vez que repitamos el experimento obtendremos un resultado diferente. De la misma forma cree que la hipótesis será cierta o falsa para el proceso estudiado pero que debido a la naturaleza aleatoria de los datos nuestro análisis puede señalarla como falsa cuando es cierta (falso negativo) o verdadera cuando es falsa (falso positivo). Por ello el enfoque frecuentista nos proporcionará, además de una respuesta a la pregunta de si la hipótesis es cierta o falsa, una probabilidad denominada p-valor. Este número es un indicador cuantas veces tendríamos que repetir la toma de datos y el análisis estadístico para obtener el resultado opuesto. Nótese que, al ser los datos aleatorios, siempre puede ser que acabe usando un conjunto que me apunte en la dirección contraria a la que realmente se comporta la población total.

Bayesiana: Hace un anális basado en la certeza que se tiene de que un hecho ocurrirá. Las observaciones se usan para actualizar o inferir la posibilidad de una hipótesis. Es una estrategia complementaria a la frecuentista, ya que con cada iteración añadiríamos un conocimiento obtenido en esa iteración para enriquecer el cálculo de la probabilidad.

USAREMOS ESTE MÉTODO CUANDO SEPAMOS A PRIORI LAS PROBABILIDADES DE CADA SUCESO, POR EJEMPLO AL LANZAR UNA MONEDA PERFECTA, POR ESO SE LLAMA CONOCIMIENTO A PRIORI,SABEMOS ANTES QUE LA PROBABILIDAD ES 50%

Por el contrario, la interpretación bayesiana se basa en un conocimiento limitado de las cosas. Afirma que sólo asocias una probabilidad a un evento porque hay incertidumbre sobre el mismo, es decir, porque no conoces todos los hechos. En realidad, un evento dado, o bien ocurrirá (probabilidad=100%) o bien no ocurrirá (probabilidad=0%). Cualquier otra cosa es una aproximación que hacemos del problema a partir de nuestro conocimiento incompleto del mismo. El enfoque bayesiano se basa por lo tanto en la idea de refinar predicciones a partir de nuevas evidencias. Un bayesiano definiría probabilidad como la expresión matemática que mide el nivel de conocimiento que tenemos para hacer una predicción. Por lo tanto, para un bayesiano, estríctamente hablando es incorrecto decir predigo que hay un 30% de probabilidades de que ocurra el evento P, sino que debería decir basándome en el conocimiento actual que tengo, tengo un 30% de certeza de que P ocurrirá.

Los frecuentistas piensan en los parámetros como unos valores fijos pero desconocidos. La estimación se basa en la elección de aquellos valores de los parámetros que maximizan la probabilidad de observar los datos.

Por su parte, los bayesianos interpretan los parámetros como variables aleatorias cuya distribución de probabilidad es estudiada en base al Teorema de Bayes. La idea es simple, un bayesiano ha de tener una distribución subjetiva de los parámetros antes de ver los datos (a priori) que modificará en función de los datos que haya observado para obtener una distribución a posteriori que resumirá todo el conocimiento del investigador sobre los parámetros de interés dados los datos y sus creencias a priori.

Los bayesianos le dan la vuelta al enfoque frecuentista. Para ellos los datos son fijos, no aleatorios porque ¿no se han obtenido ya? ¿Qué tienen de aleatorios los valores anotados en esas tablas que hemos ido generando? En cambio la hipótesis es para ellos aleatoria y puede ser o no verdad. Los bayesianos no determinan si la hipótesis es cierta o falsa, únicamente calculan la probabilidad (en sentido bayesiano) de que sea cierta o falsa.

Finalmente, sabemos que hay muchos resultados frecuentistas que pueden obtenerse desde una perspectiva bayesiana aunque la filosofía subyacente sea diferente. Mínimos cuadrados ordinarios (MCO) es el ejemplo más sencillo. Dicho estimador frecuentista coincide exactamente con la media de la distribución bayesiana bajo unas creencias concretas en el marco del modelo lineal normal. Existen en la literatura muchos otros ejemplos.

El objetivo de las técnicas de inferencia estadística es contrastar la hipótesis nula (Ho) y la hipótesis alternativa (H1), y en este sentido se pueden cometer dos tipos de error:

Error tipo I o alfa --> Rechazar la hipótesis nula siendo esta verdadera.
Error tipo II o beta --> Rechazar la hipótesis alternativa siendo esta cierta. La probabilidad complementaria del tipo II se llama poder de contraste (p-value) y es 1-beta.Este valor está totalmente impactado por el tamaño de la muestra, ya que el intervalo de confianza se va reduciendo cuanto mayor es el espacio muestral, por tanto más probable es rechazar Ho de manera efectiva.

La base de la inferencia estadistica es la probabilidad. Esta, asigna un número entre 0 y 1 que significa la posibilidad de que sea cierto. La probabilidad es el modelo de inferencia que se usa habitualmente para hacer estimaciones sobre fenómenos aparentemente aleatorios.

Teniendo en cuenta la conexión entre ambos enfoques en la práctica y la posibilidad de obtener los mismos resultados bajo ambos enfoques, los bayesianos argumentan que su interpretación del problema es siempre más intuitiva y natural. Recordemos que un bayesiano proporcionaría conclusiones del tipo: hay un 95% de probabilidad de que el parámetro esté entre 0,3 y 0,8, sin embargo, un frecuentista afirmaría: si generáramos 100 muestras aleatorias del mismo tamaño y repitiéramos la estimación 100 veces, en 95 de ellas el parámetro estimado estaría entre 0,3 y 0,8.

Los contrastes de hipótesis parecen por tanto más naturales en el marco bayesiano. Para un bayesiano convencido, un p-valor o un t-ratio es una herejía, ya que lo único que necesitamos para contrastar hipótesis es tener la distribución a posteriori de los parámetros. Salvo en contadas excepciones, un frecuentista basa su contraste de hipótesis en el análisis asintótico, es decir, en calcular p-valores de la distribución asintótica del estimador (no del parámetro verdadero, que es un valor fijo). Esta distribución, mayoritariamente normal gracias a numerosos teoremas centrales del límite, es la que el estimador tendría si el investigador tuviera muchos más datos de los que realmente tiene. ¿es ésta una buena forma de hacer inferencia en general? Por el contrario, el bayesiano se basa en la distribución de los parámetros dados sus datos, que puede en principio tener cualquier forma no gaussiana.

Al admitir información adicional, los métodos bayesianos reconocen que cada problema es distinto y promueven que el procedimiento de análisis se adapte al problema en cuestión, y no al revés; como consecuencia, tienden a ser más flexibles. Sin embargo, para poder hacer esto posible, la implementación de las técnicas bayesianas usualmente requiere de un esfuerzo computacional muy alto. La mayor parte de este esfuerzo se concentra en el cálculo de ciertas características de la distribución final del parámetro de interés. Por ejemplo, hay que integrar para pasar de una distribución conjunta a una colección de distribuciones marginales que sean útiles para hacer inferencias sobre los parámetros de interés. En la gran mayoría de los problemas las integrales requeridas no pueden resolverse analíticamente, por lo que es necesario contar con métodos numéricos eficientes que permitan calcular o aproximar integrales en varias dimensiones.

El enfoque bayesiano no sustituye el estudio de las frecuencias, lo enriquece con experiencias anteriores ya que puede ir refinando hacia atrás los resultados obtenidos en cada iteración. Un ejemplo extremo es que si en los últimos 4.500 millones de años en Cuenca ha salido el sol por el levante, suponemos que la probabilidad de que mañana vaya a volver a salir el sol es más alta con un estudio bayesiano de lo que sería con un análisis de estadística clásica, que no da ningún peso a la historia pasada. Pero cuando un médico pronostica le quedan unos seis meses de vida, querríamos creer que el galeno nos está facilitando esta información no sólo en base a un análisis de las frecuencias para esta dolencia, sino que además incorpora su experiencia de otros diagnósticos similares, y en especial nuestro propio historial clínico. Las conclusiones son afectadas por la premisa, pero si ésta es bien formulada, más que un defecto sería una ventaja

Afortunadamente, durante la segunda mitad del siglo XX se comenzaron a desarrollar técnicas numéricas flexibles y eficientes basadas en métodos de simulación estocástica. Esto, aunado al desarrollo de la tecnología que ha dado lugar a una mayor capacidad de procesamiento y de almacenamiento de los equipos de cómputo gracias a la virtualización, el Cloud Computing y el Big Data, está provocando un auge de los métodos bayesianos.

Probablidad es la cantidad de población que resume una aleatoriedad.Dos herramientas muy útiles para el cálculo de probabilidad son la distribución de probabilidad y la función de distribución

Esto es, ante la pregunta de si una hipótesis es correcta, los frecuentistas responden un sí o un no pero apostillan con qué frecuencia obtendrías el resultado opuesto mientras que los bayesianos no responden directamente y se limitan a indicar cómo de seguros están de uno u otro resultado.

Haremos usos de ambos métodos dependiendo de lo que necesitemos calcular. Usaremos por ejemplos métodos frecuentistas para el cálculo del error muestral y comparación entre muestras. Mientras que usaremos métodos bayesianos para el cálculo de especificidad y sensibilidad.

TEORÍA DE LA PROBABILIDAD

Dado un experimento aleatorio (por ejemplo lanzar un dado) su medida de probabilidad es la cantidad de dicha población que mide la aleatoriedad.

Está muy bien explicado en el siguiente artículo, del que he tomado algunas notas:

https://blogs.20minutos.es/mati-una-profesora-muy-particular/tag/probabilidad-condicionada/

Específicamente, la probabilidad toma un posible valor del experimento dónde:

• Assigna un número entre 0 y 1.
• La probabilidad de que algo ocurra seguro es 1.
• La probabilidad de la unión de dos conjuntos de posibles resultados que no tienen nada en común (son mutuamente excluyentes, por ejemplo que salga un 2 o un 4 en un lanzamiento) es la suma de sus respectivas probabilidades.

Para asignar una probabilidad a un suceso, echamos mano de la regla de Laplace

laplace

En un caso más ámplio, como en el lanzamiento de dados

laplace

El matemático ruso Kolgomorov formalizó estas reglas de probabilidad.

• La probabilidad de nada ocurra es 0
• La probailidad de que sí ocurra es 1
• La probabilidad de un hecho es 1 menos la probabilidad del hecho opuesto

kolgomorov

• La probabilidad de que sucedan al menos uno de dos (o más) cosas que no pueden ocurrir simultáneamente es la suma de sus probabilidades (UNIÓN DE CONJUNTOS)

kolgomorov

• Si un evento A implica la ocurrencia de un evento B, entonces la probabilidad de que ocurra A es menor que la probabilidad de que ocurra B (INCLUSIÓN DE CONJUNTOS)

qownnotes-media-GBLpWb

• Para cualquiera de dos eventos (o más) la probabilidad de que ocurran es la suma de sus probabilidades menos su intersección. (INTERSECCIÓN DE CONJUNTOS)

qownnotes-media-waZsrq

kolgomorov

qownnotes-media-mTAmEg

\[ \begin{eqnarray*} A_1 & = & \{\mbox{Person has sleep apnea}\} \\ A_2 & = & \{\mbox{Person has RLS}\} \end{eqnarray*} \] \[ \begin{eqnarray*} P(A_1 \cup A_2 ) & = & P(A_1) + P(A_2) - P(A_1 \cap A_2) \\ & = & 0.13 - \mbox{Probability of having both} \end{eqnarray*} \]

Una Variable aleatoria es el resultado numérico de un experimento. Las variables aleatorias tienen dos naturalezas, o discreta o contínua.

Un ejemplo de variable discreta es todo lo que sea contar un número de cosas(por ejemplo machos y hembras, o cualquier otra variable categórica), un ejemplo de continua es la medida de tu peso durante un mes.

Supongamos que deseamos conocer algún parámetro poblacional, por ejemplo la edad media de una población. En el contexto de la estadística clásica (frecuentista), podríamos extraer al azar múltiples y diferentes muestras representativas de forma sucesiva y determinar el estimador más apropiado en cada una de esas muestras. Para este caso, el estimador más apropiado para la media poblacional es la media muestral. Cada una de las medias de estas muestras tomarían un determinado valor diferente para cada muestra (o lo que es lo mismo, el estimador es tratado como una variable en sí mismo).En estadística, un estimador es un estadístico (esto es, una función de la muestra) usado para estimar un parámetro desconocido de la población.

Suponiendo la normalidad de la variable edad, los valores que irían tomando estas medias muestrales se distribuirían según una curva de Gauss (algunos de estos valores se darían con mayor frecuencia y otros con menor frecuencia). El área bajo dicha curva de distribución de frecuencias del estimador es la probabilidad total, es decir, la unidad.

DISTRIBUCIÓN DE PROBABILIDAD

Función de probabilidad (función de masa de probabilidad, PMF en ingles DISCRETOS)

Función que determina la probabilidad de cada evento posible, cuando los posibles valores son DISCRETOS: ejemplo:

• valor p(1) -> 0.3
• valor p(4) ->0.6
• valor p(7)->0.1

PMF valores discretos

Las distribuciones de probabilidad en base a su PMF más comunes son Bernoulli (valor p o p-1 muy útil para la prevalencia de una enfermedad), binomial (bernuilli repetido n veces) o poisson(contar el número de veces que se da un evento en el tiempo)

PMF valores discretos

Función de densidad (PDF en ingles CONTINUO)

Función que determina la probabilidad de cada evento posible, cuando los posibles valores son CONTINUOS, por ejemplo probabilidad asociada a tiempo de espera en urgencias.

Ejemplo la proporción de llamadas que pueden ser atendidas en un día para dar citas en consultas externas sigue la siguiente curva:

\[ f(x) = \left\{\begin{array}{ll} 2 x & \mbox{ for } 1 > x > 0 \\ 0 & \mbox{ otherwise} \end{array} \right. \]

x <- c(-0.5, 0, 1, 1, 1.5); y <- c( 0, 0, 2, 0, 0)
plot(x, y, lwd = 3, frame = FALSE, type = "l")

Cuál será la probabilidad de atender el 75% al menos de las llamadas? Esta gráfica representa su PDF

Esta operación es muy sencilla ya que se puede calcular con trigonometría o através de la función beta, que es justo para distribuciones de probabiidad trianangular

1.5 * .75 / 2

## [1] 0.5625

pbeta(.75, 2, 1)

## [1] 0.5625

Por tanto la probabilidad de atender al menos el 75% de las llamadas es del 56%

Ya que calculamos el área bajo la curva, si elegimos un valor concreto su pdf nos dara 0,ya que el area de una recta es 0.

Area bajo curva de una recta es 0

Las distribuciones de probabilidad en base a su PDF más comunes son Normal (la normal estandar es la que tiene e(x) = 0 y d(x) = 1)o de Gaus o exponencial (el tiempo entre dos eventos)

Los histogramas, pueden representar esta función de densidad en base a probabilidades o bien en base a número de ocurrencias de cada evento. EN nuestro caso, esto es fundamental ya que las probabilidades que vamos a manejar en el mundo de la la sanidad están todas basadas en la ocurrencia de un suceso, no en un cálculo a priori de sus probabilidades.

A la hora de representarlo, el histograma será mucho menos fiable para conocer su PDF ya que depende en gran medida del número de separaciones que hagamos, es más fiable utilizar diagrama de densidad kernel

http://www.statmethods.net/graphs/density.html

pbeta(c(0.4, 0.5, 0.6), 2, 1)

## [1] 0.16 0.25 0.36

sqrt(0.5)

## [1] 0.7071068

qbeta(0.5, 2, 1)

## [1] 0.7071068

PROBABILIDAD CONDICIONAL

La probabilidad de que un suceso pase sabiendo que ha pasado otro. En este caso, el hecho de conocer de antemano un hecho relaccionado (DEPENDIENTE) nos hace modificar la probabilidad de que ocurra el hecho en el que estamos interesados.

Esta es la base de la aproximación Bayesiana, ya que aportamos conocimiento a nuestra predicción y modificamos su probabilidad respecto a ese conocimiento.

• Pongamos un hecho B = de manera que \(P(B) > 0\)
• La probabilidad condicionada de otro evento A, sabiendo previamente que ha ocurrido B es definido como la probabilidad de que los DOS ocurran dividido por la probabilidad de que B occurra \[P(A\:|\:B) = \frac{P(A \:\cap\: B)}{P(B)}\]
• En el caso de A y B sean independientes, entonces \[P(A\:|\:B) = \frac{P(A)P(B)}{P(B)} = P(A)\]

Un ejemplo de probabilidad condicionada, es qué probabilidades tengo de sacar un 1 si sé que ha salido impar.

qownnotes-media-jvqfcq

El típico ejemplo de sucesos dependientes es el de sacar bolas de una bolsa. Si tenemos 5 bolas rojas y 3 bolas verdes en una bolsa, ¿cuál es la probabilidad de meter la mano y sin mirar sacar una bola roja?

sucesosDependientes

La dejamos fuera ¿Y la probabilidad de que la segunda sea roja?Esto implica que el resultado del segundo experimento depende del primero. Lo que queremos conocer es la probabilidad de que en la segunda extracción (ya hay una bola fuera) saquemos una bola roja. A esa probabilidad que queremos calcular la llamaremos P(2ªR). Para calcular P(2ªR) vamos a necesitar calcular otras probabilidades antes. Por ejemplo, la probabilidad de que la 2ª bola sea roja condicionada al hecho de que la primera que salió fuese roja, lo llamamos P(2ªR | 1ªR)

sucesosDependientes

Una vez tenemos todas las partes implicadas, será sumar las probabilidades de que la segunda haya sido roja condicionada a que la primera sea roja y que la segunda sea roja condicionando a que la primera no sea roja

Este caso es muy sencillo y podemos haer los cálculos símplemente contando

sucesosDependientes

Ya sólo queda hacer la suma:

sucesosDependientes

Teorema de Bayes Es el uso más famoso de la probabilidad condicional. Calcula la probabilidad de A condicionado a B desde la probabilidad de que pase B condicionado a A. Este teorema nos da la posibilidad de conocer las probabilidades condicionadas cruzadas de ambos eventos entre sí.

\[P(B\:|\:A) = \frac{P(A\:|\:B)P(B)}{P(A\:|\:B)P(B)+P(A\:|\:B^c)P(B^c)}\] donde \(B^c\) = corresponde a la probabilidad de que no suceda el evento \(B\), \(P(B^c) = 1 - P(B)\)

Es muy útil para calcular la efectividad de test diagnósticos ya que nos relacciona la efectividad del resultado de un test (Hecho A) con que realmente se tenga la enfermedad (Hecho B). Ejemplo: + y - son los resultados de diagnóstico de un test médico, y D que el sujeto tenga la enfermedad.

Probabilidad de que el test sea positivo = \(P(+)\) Probabilidad de que el test sea negativo = \(P(-)\) Probabilidad de que se tenga la enfermedad = \(P(D)\) Probabilidad de que no se tenga la enfermedad = \(P(D^c)\)

SENSIBILIDAD es la probabilidad de que el test sea positivo sabiendo que individuo tiene la enfermedad (Nos indica cómo ha funcionado). Es desable que sea alto, ya que indicará que nuestro test es fiable cuando da positivo.

Sensibilidad = \(P(+\:|\:D)\)

ESPECIFIDAD es la probabilidad de que el test sea negativo sabiendo que el individuo no tiene la enfermedad (Nos indica cómo ha funcionado). También es desable que sea alto, ya que indicará que nuestro test es fiable cuando da negativo.

Especificidad = \(P(-\:|\:D^c)\)

Prevalencia aparente Es la probabilidad de tener sepsis observada = \(P(D)\)

Prevalencia real es practicamente imposible de determinar pero se puede aproximar desde la prevalencia aparente con el estimador Rogan-Gladen

Valor predictivo positivo Probabilidad de que tenga la enfermedad sabiendo que el test es positivo (Nos indica cómo va a funcionar, tiene en cuenta prevalencia) = \(P(D\:|\:+)\)

Valor predictivo negativo Probabilidad de no tener la enfermedad sabiendo que el test es negativo (Nos indica cómo va a funcionar, tiene en cuenta prevalencia) = \(P(D^c\:|\:-)\)

Razón de verosimilitud diagnóstica -> Existen con el fin de dimensionar el beneficio clínico de un test de manera independiente a la prevalencia, ya que esta puede variar.

Razón de verosimilitud diagnóstica de test positivo El ratio entre la posibilidad de que en un paciente séptico haya saltado la alerta y que en un paciente no séptico haya saltado la alerta \[DLR_+ = \frac{sensitivity}{1-specificity} = \frac{P(+\:|\:D)}{P(+\:|\:D^c)}\] \(DLR_+ = .997 / (1 - .985) \approx 66\) DLR+ Ratio de 66 significa que la probabilidad de tener la enfermedad después del test es 66 veces mayor que antes del test

Razón de verosimilitud diagnóstica de test negativo Es el ratio entre la posibilidad de que en un paciente séptico no haya saltado la alerta y la posibilidad de que en un paciente no séptico no haya saltado la alerta. \[DLR_- = \frac{1 - sensitivity}{specificity} = \frac{P(-\:|\:D)}{P(-\:|\:D^c)}\] \(DLR_- = (1 - .997) / .985 \approx .003\) DLR- es 0.005 significa que la probabilidad de no tener la enfermedad antes del test es 0.003 veces mayor que después del test

Odds ratio Asociación entre exposición y resultado. Ejemplo de razón de oportunidad de indice de cancer de pulmón en fumador es 2. Esto significa que fumadores tienen 2 veces más probab. de tener cancer. = \(p/(1-p)\)

Diferencia entre probability y odds:

Probabilidad 0.80 Odds 0.80/1-0.80 = 4 –> Odds de 4 a 1

* $\frac{P(D)}{P(D^c)}$ = **pre-test odds**, o razón de oportunidad de enfermedad en ausencia del test
* $\frac{P(D\:|\:+)}{P(+\:|\:D^c)}$ = **post-test odds**, o razón de oportunidad de enfermedad dado un test positivo
* $DLR_+$ = factor por el cual la razón de oportunidad en la presencia de un test positivo puede ser multiplicado para obtener el post-test odds
* $DLR_-$ = relacciona la disminución en la razón de oportunidad de enfermedad después de un test negativo
* Siguiendo el ejemplo, para una sensibilidad de 0.997 y especificidad de 0.985, los razones de verosimilitud diagnosticos son $$DLR_+ = .997/(1-.985) = 66 ~~~~~~ DLR_- =(1-.997)/.985 = 0.003$$
* Esto indica que en el resultado de un test positivos, la razón de verosimilitud de la enfermedad es 66 veces el pretest odds

Prevalencia aparente= (8503+37)/(8503+37+2758+355024)=0.023

Prevalencia real = (Prev.apa+Espec-1)/(Sensib.+Espec-1) = (0.023+0.992-1)/(0.995+0.992-1)= 0.01519757

Proporción falso positivo = 2758/11261 = 0.24

Proporción falson negativo = 37/355061 = 0.0001

Sensibilidad = 8503/(8503+37) = 0.995

Especificidad = 355024/(2758+355024) = 0.992

VPP = 8503/(8503+2758) = 0.75

VPN = 355024/(2578+355024) = 0.992

DLR+ = (sensibilidad/1-especificidad) = 124

DLR- = (1-sensibilidad/especificidad) = 0,005

Odds ratio (razón de oportunidades) = (8503/2758)/(37/355024) = 29582.43

Ratio o razón –> Indice obtenido al dividir dos cantidades. Ninguno de los elementos del numerador estón incluidos en el denominador

Ejemplo Ratio alerta post = 8503/2758 = 3.08 por cada 3 alertas salta 1 que no es

Proporción -> El numerador está incluido en el denominador

    Proporcion alerta post =  8503/(8503+2758) = 0.75

Tasa -> Magnitud de cambio de un parámetro por unidad de cambio de otro que incluye una medida adiconal de tiempo en el denominador

Tasa de pacientes septicos en 2018: 8540/366322 = 0.023 es una tasa
de 23 pacientes septicos por cada 1000 urgencias.

qownnotes-media-tBAWRb

POSITIVE PREDICTIVE VALUE es la probabilidad de tener la enfermedad sabiedo que el test es positivo.

NEGATIVE PREDICTIVE VALUE es la probabilidad de no tener la enfermedad sabiendo que el test es negativo.

En la ausencia de un test, se llama PREVALENCIA a la probabilidad de tener esa enfermedad.

qownnotes-media-dOUpkP

qownnotes-media-tLsDNl

qownnotes-media-pPYGcB

qownnotes-media-RUIsxz

DLR+ es cada cuantos test positivos acertados se da un test negativo err?neo en este caso 998/15 = 66, cada 66 test positivos acertados se da 1 negativo err?neo

DRL- es cada cuantos test positivos erróneos se da un test negativo acertado.

qownnotes-media-PaNiRC

qownnotes-media-xTKrgI

qownnotes-media-CtTfPs

qownnotes-media-DynpEw

qownnotes-media-htWhsk

qownnotes-media-EtEyvi

ESTUDIO DETALLADO DE LA SENSIBILIDAD Y ESPECIFICIDAD

Una vez que la prueba se ha aplicado a una población determinada, varios índices se utilizan para evaluar la prueba: - Verdadero positivos (VP): Número de casos que la prueba declara positivos y que son verdaderamente positivos. - Falsos positivos (FP): Número de casos que la prueba declara positivos y que en realidad son negativos. - Verdadero negativo (VN): Número de casos que la prueba declara negativos y que son realmente negativos. - Falsos negativos (FN): Número de casos que la prueba declara negativos y que en realidad son positivos. - N = TP + FP + FN + TN

qownnotes-media-HqjvcR

• Sensibilidad (equivalente a la tasa de positivos verdaderos): Proporción de casos positivos que están bien detectadas por la prueba. La definición matemática es: Sensibilidad = VP / (VP + FN).

• Especificidad (también llamada Tasa de verdaderos negativos): proporción de casos negativos que son bien detectadas por la prueba. La definición matemática es: Especificidad = VN / (VN + FP).

• Tasa de falsos positivos (FPR): Proporción de casos negativos que la prueba detecta como positivos.

• Falso Negativo Precio (FNR): Proporción de casos positivos que la prueba detecta como negativo. - Prevalencia: la frecuencia relativa de los acontecimiento de inter?s en la muestra total (VP + FN) / N.

• Valor Predictivo Positivo (VPP): Proporción de casos verdaderamente positivos entre los casos positivos detectados por la prueba. Tenemos PPV = TP / (TP + FP), o PPV = Sensibilidad x Prevalencia / [(Sensibilidad x Prevalencia + (1-Especificidad) (1-Prevalencia)]. Es un valor fundamental que depende de la prevalencia, un índice que es independiente de la calidad de la prueba.

• Valor predictivo negativo (VPN): Proporción de casos verdaderamente negativos entre los casos negativos detectados por la prueba. Tenemos VPN = VN / (VN + FN), o VPP = Especificidad x (1 - Prevalencia) / [(Especificidad (1-Prevalencia) + (1-sensibilidad) x Prevalencia] . Este ?ndice depende tambi?n de la prevalencia que se independiente de la calidad de la prueba.

• Razón de verosimilitud positiva (LR +): Esta relación indica a qué punto una persona tiene más posibilidades de ser positivo en la realidad cuando la prueba es que está diciendo es positivo. Tenemos LR + = sensibilidad / (1-especificidad). La RP + es un valor positivo o nulo.

• Razón de verosimilitud Negativa (LR-): Esta relación indica a qué punto una persona tiene más posibilidades de ser negativo, en realidad, cuando la prueba es que está diciendo es positivo. Tenemos LR-= (1-sensibilidad) / (especificidad). El LR-es un valor positivo o nulo.

• Odds ratio: El odds ratio indica la cantidad de un individuo es más probable que sea positivo si el resultado es positivo, en comparación con los casos en que la prueba es negativa. Por ejemplo, un odds ratio de 2 significa que la probabilidad de que el caso positivo se produce es dos veces superior si la prueba es positiva que si es negativo. El odds ratio es un valor positivo o nulo. Tenemos Odds ratio = VPxVN / (FPxFN).

• El riesgo relativo: El riesgo relativo es un ratio que mide mejor la prueba se comporta cuando se trata de un informe positivo que cuando es negativo. Por ejemplo, un riesgo relativo de 2 significa que la prueba es dos veces más potente cuando es positiva que cuando es negativo. Un valor cercano a 1 corresponde a un caso de independencia entre las filas y columnas, y una prueba de que funciona tan bien cuando es positiva que cuando es negativo. El riesgo relativo es un valor nulo o un valor positivo dado por: riesgo relativo = VP / (VP + FP) / (FN / (FN + VN)).

VARIABLES INDEPENDIENTES

En este apartado cambiamos totalmente de premisa, ya que en los anteriores apartados las variables estaban relaccionadas entre sí, en este caso que ocurra la primera no infiere en que ocurra la segunda.

Dos variables son independientes si la ocurrencia de una no tiene nada que ver en la ocurrencia de otra.

Por ejemplo tirar dos veces una moneda y que salga las dos cara, son sucesos independientes 0.5x0.5 = 0.25

Si los sucesos no son independientes la intersección de las probabilidades no se pueden calcular multiplicando.

• dos eventos \(A\) y \(B\) son independientes si lo siguiente se cumple:
  • \(P(A\:\cap\:B) = P(A)P(B)\)
  • De manera equivalente \(P(A\:|\:B) = P(A)\)
• dos variables \(X\) e \(Y\) son independientes, si para cualquier subconjunto, A y B, se cumple
  • \(P([X \in A]\cap[Y \in B]) = P(X \in A)P(Y \in B)\)
• independence = statistically unrelated from one another
• si \(A\) es independiente de \(B\), entonces se cumple
  • \(A^c\) es independiente de \(B\)
  • \(A\) es independiente de \(B^c\)
  • \(A^c\) es independiente de \(B^c\)

Son sucesos independientes, por ejemplo, el lanzamiento consecutivo de monedas. Si empezamos en el punto azul y caminamos sobre las aristas azules, tenemos todos los resultados

SucesosIndependientes

Si queremos saber qué probabilidad hay de sacar 4 caras seguidas, no tenemos más que seguir el camino

SucesosIndependientes

Se trata de calcular la probabilidad de una intersección de sucesos independientes, porque el resultado de cada lanzamiento es independiente del anterior. Un suceso es que la primera salga cara, otro que la segunda salga cara, el tercero que salga cara en la tercera y el cuarto, que la cuarta tirada dé cara. Cada uno de estos sucesos, tiene probabilidad 1/2, como hemos dicho, Eso significa que en cada tirada de las 4, la probabilidad de que salga cara es 1/2 o sea, de un 50%.

SucesosIndependientes

Ahora bien, la probabilidad de que las cuatro veces nos salga cara es la probabilidad de la intersección, y eso se calcula sabiendo que, como son independientes, la probabilidad de la intersección es igual al producto de las probabilidades.

SucesosIndependientes

Es muy fácil caer aquí en la falacia del apostador o de Montecarlo. Este nombre viene porque en su casino, sucedió un hecho insólito hasta la fecha, una noche en la ruleta salieron 20 números negros consecutivos, por lo que las apuestas por roja en la siguiente tirada fueron en masa, con la creencia de que si habían salido 20 resultados seguidos negros, el siguiente sería muy difícil que fuera negro. Esa noche llegaron a salir hasta 26 veces seguidas números negros. Es la falsa creencia de que lo sucedido anteriormente en un experimento aleatorio, al azar, afectará al resultado de los experimentos siguientes.

Se están confundiendo la probabilidad de obtener 5 caras seguidas con la probabilidad de obtener cara en el 5º lanzamiento.

SucesosIndependientes

Todos los estudios que vamos a realizar a continuación, no sólo están basados en sucesos independientes, si no que además están identicamente distribuidos, IID variables

• se dice que dos variables aleatorias son IID si son independientes e identicamente distribuidas
  • independientes = estadisticamente no relaccionadas entre sí
  • identicamente distribuidas = ambas siguen la misma distribución poblacional.
• variables aleatorias IID = modelo por defecto para muestras aleatorias = punto inicial por defecto para el estudio de inferencia

VALOR ESPERADO O ESPERANZA MATEMATICA(Centro geométrico de las Áreas de probabilidad) o MEDIA

Se insiste mucho en la diferencia entre muestra y población. En este caso cuando hagamos cálculos sobre una muestra haremos los cálculos en base a los valores encontrados en esa muesta (Ej:tiro un dado 10 veces y salen 2 1s, 3 2s, 1 3, 4 5s.). Cuando se habla sobre los cáculos sobre una población estaremos hablando sobre las probabilidades que tenga una variable de tomar un valor sin haber realizado muestras ( o lo que es lo mismo si hubiéramos hecho infinitas muestras)

MEDIA DE UNA POBLACIÓN

Es una propiedad de las distribuciones de probabilidad de la misma manera que lo son la mediana (valor del punto que está en la mitad) y la varianza

En variables discretas será el valor de su probabilidad multiplicado por el propio valor

E(x)= x1P(x1)+x2P(x2)+…

\[E[X] = 1 \times \frac{1}{6} + 2 \times \frac{1}{6} + 3 \times \frac{1}{6} + 4 \times \frac{1}{6} + 5 \times \frac{1}{6} + 6 \times \frac{1}{6} = 3.5\]

qownnotes-media-S12388

En variables contínuas, tenemos que hacer uso de la función de distribución para recoger la probabilidad acumulada.(dx)

E(x) = Integral (xf(x) d(x))

En PDF es lo mismo , el valor esperado es el centro de masa del área que encierra su función de densidad. Una propieda de la media es que su valor es el mismo en la función de densidad (DE POBLACION) que en la función de densidad de la media de X muestras del mismo tipo (MUESTRA) , es decir, si el valor esperado de una distribución normal centrada en 0 es 0 (en azul), función de densidad de la media de 10 distribuciones normales distintas será también 0 (en rosa). Cuanto mayor sea la población sobre la que estamos sacando la media (en lugar de 10 que fuera 100) más se aglutinarán las densidades entorno al valor esperado.

qownnotes-media-z14800

VARIABILIDAD

VARIANZA DE UNA MUESTRA

Es la suma de la distancia de todos los puntos de la muestra respecto a la media (se pone al cuadrado para evitar distancias negativas y calcular su distancia absoluta)

S^2 = (1-3.9 +1-3.9 +5-3.9 +4-3.9 +3-3.9 +4-3.9 +6-3.9 +5-3.9 +4-3.9 +6-3.9 )^2 / 10 -1 = 0

S = Desviación estándar Se aplica a raiz cuadrada para que esté en las mismas unidades que el resto.

La varianza de una muestra se define como \[S^2 = \frac{\sum_{i=1} (X_i - \bar X)^2}{n-1}\]

También es una variable aleatoria:

• Tiene asociada una distribución poblacional
• Su valor experado es la varianza poblacional
• Su distribución se concentra cada vez más alrededor de la varianza poblacional cuantos más datos haya.
• Su raíz cuadrada es la desviación estándar

En el caso de la varianza de las muestras hay una pequeña variación ya que se divide por n-1. Esto es debido a que las primeras n-1 observaciones son son independientes de la media. la úiltima no es independiente ya que puede ser calculada por medio de la media de la muestra usada en la formula. En otras palabras, la varianza de la muestra es CASI la media de la desviación al cuadrado de la media de la muestra.

Var(x) suma var(x) /n-1

La distribución de la varianza de las muestras también se verá afectado por el tamaño de la muestra, de manera que a mayor tamaño de la muestra la curva se irá acercando a la variaza de la población que está intentando estimar:

qownnotes-media-M19368

Var(X’)=Var(1/n*Sum(X_i))=(1/n^{2)Var(Sum(X_i))=(1/n^2)Sum(sigma}2)=sigma^2/n para poblaciones infinitas cuando las muestras son independientes | infinite populations when our samples are independent

MEDIA Y DESVIACIÓN ESTANDAR (ERROR ESTANDAR) DE MÚLTIPLES MUESTRAS EN DISTINTAS DISTRIBUCIONES

Si simulamos con R 1000 muestras de una distribción normal con 10 valores posibles vemos que se acerca mucho a lo calculado para la población, y cuanto más se aumenten las muestras más cerca estará ese valor:

# specify number of simulations with 10 as number of observations per sample
nosim <- 1000; n <-10
# estimated standard deviation of mean
sd(apply(matrix(rnorm(nosim * n), nosim), 1, mean))

## [1] 0.3037642

# actual standard deviation of mean of standard normals
1 / sqrt(n)

## [1] 0.3162278

# estimated standard deviation of the sample means
sd(apply(matrix(runif(nosim * n), nosim), 1, mean))

## [1] 0.09341847

# actual standard deviation of the means
1/sqrt(12*n)

## [1] 0.09128709

# estimated standard deviation of the sample means
sd(apply(matrix(rpois(nosim * n, lambda=4), nosim), 1, mean))

## [1] 0.6606151

# actual standard deviation of the means
2/sqrt(n)

## [1] 0.6324555

# estimated standard deviation of the sample means
sd(apply(matrix(sample(0 : 1, nosim * n, replace = TRUE), nosim), 1, mean))

## [1] 0.1601838

# actual standard deviation of the means
1/(2*sqrt(n))

## [1] 0.1581139

EJEMPLOS CON DISTRIBUCIONES

BERNULLI

Si X es una variable aleatoria que mide el número de éxitos, y se realiza un único experimento con dos posibles resultados (éxito o fracaso), se dice que la variable aleatoria X se distribuye como una Bernoulli de parámetro p. (Ejemplo tirar 1 moneda 1 vez y que salga cruz)

• Distribución de Bernoulli es el resultado de un experimento con sólo 2 resultados (binario) * 1 = “éxito” con probabilidad de \(p\) * 0 = “fallo” con probabilidad de \(1 - p\)
  • PMF se define como \[P(X=x) = p^x(1 - p)^{1-x}\]
  • mean = \(p\)
  • variance = \(p(1 - p)\)

Si hacemos varias observaciones de Bernuilli pongamos \(x_1,\ldots, x_n\) , la posibilidad (likelihood)de que el evento sea éxito es \[ \prod_{i=1}^n p^{x_i} (1 - p)^{1 - x_i} = p^{\sum x_i} (1 - p)^{n - \sum x_i} \] \(\hat p = \sum_i x_i / n\) es el máximo estimador de posibilidad para \(p\)

Todas las posibilidades para n pequeña

BINOMIAL

cuenta el número de éxitos en una secuencia de n ensayos de Bernoulli independientes entre sí, con una probabilidad fija p de ocurrencia del éxito entre los ensayos. Un experimento de Bernoulli se caracteriza por ser dicotómico, esto es, solo dos resultados son posibles

Para n = 1, la binomial se convierte, de hecho, en una distribución de Bernoulli.

• variable aleatoria binomial = suma de n variables de Bernoulli \[X = \sum_{i=1}^n X_i\] donde \(X_1,\ldots,X_n = Bernoulli(p)\)
  • PMF se define como \[P(X=x) = {n \choose x}p^x(1-p)^{n-x}\] donde \({n \choose x}\) = número de maneras de seleccionar \(x\) veces de \(n\) opciones sin reemplazamiento \(x=0,\ldots,n\)
  • combinación o “\(n\) combinaciones de \(x\)” se define como \[{n \choose x} = \frac{n!}{x!(n-x)!}\]
  • Los casos base son \[{n \choose n} = {n \choose 0} = 1\]

* de 8 hijos, cuál es la probabilidad de tener 7 chicas o más (50/50 chance)? \[{8 \choose 7}.5^7(1-.5)^{1} + {8 \choose 8}.5^8(1-.5)^{0} \approx 0.04\]

• choose(8, 7) = función de R function para calcular las combinaciones de x en n elementos
• pbinom(6, size=8, prob =0.5, lower.tail=TRUE) = probabilidad de 6 o menos éxitos de 8 muestras con probabilidad 0.5 (CMF)
  • lower.tail=FALSE = devuelve el complemento, en este caso es la probabilidad de que sea mayor a 6 de 8 muestrasy of 0.5

qownnotes-media-fjsoAv

plot(pvals, dbinom(7, 8, pvals) / dbinom(7, 8, 7/8) ,
     lwd = 3, frame = FALSE, type = "l", xlab = "p", ylab = "likelihood")

DISTRIBUCION NORMAL ESTANDAR

Es una distribución normal pero con la particularidad de que tiene mu = 0 y sigma = 1

qownnotes-media-BjqEKS

De manera que 1 desviación estandar recoge el 68% de los datos, 2 desviaciones estándar el 95% y 3 desviaciones estandar el 99%.

qownnotes-media-SIFoSu

qownnotes-media-hjbGgE

# specify number of simulations with 10 as number of observations per sample
nosim <- 1000; n <-10
# estimated standard deviation of mean
sd(apply(matrix(rnorm(nosim * n), nosim), 1, mean))

## [1] 0.3166805

# actual standard deviation of mean of standard normals
1 / sqrt(n)

## [1] 0.3162278

DISTRIBUCION NORMAL

Por tanto una distribucion normal tendrá mu != 0 y sigma != 1, la única diferencia con los gráficos anteriores es que no estará centrada en 0, y que los ejes de las x cada 1,2 y 3 estar? multiplicado por su sigma.

DistribuciónNormal

qownnotes-media-mAfrpi

• Approximadamente en 68%,95% y 99% de la densidad normal recae en 1,2 y 3 desviaciones estándar de la media.

• -1,28,-1.645,-1.96 and -2.33 son los percentiles 10º,5º, 2.5º y 1º de la distribución

• 1,28,1.645,1.96 and 2.33 son por simetría los percentiles 90º,95º, 97.5º y 99º de la distribución

qownnotes-media-RxbLpG

qownnotes-media-R14800

qownnotes-media-B14800

qownnotes-media-P19368

qownnotes-media-L19368

POISSON

expresa, a partir de una frecuencia de ocurrencia media, la probabilidad de que ocurra un determinado número de eventos durante cierto periodo de tiempo. Concretamente, se especializa en la probabilidad de ocurrencia de sucesos con probabilidades muy pequeñas, o sucesos “raros”.

Tanto el valor esperado como la varianza de una variable aleatoria con distribución de Poisson son iguales

qownnotes-media-mMWyDG

qownnotes-media-BmoCZj

qownnotes-media-Y14800

qownnotes-media-i14800

qownnotes-media-i19368

qownnotes-media-u19368

Es muy útil para:

*Conteos de datos Cuantas veces se repite un determinado evento.

*Modelar eventos de tiempo o de supervivencia, que modeliza el tiempo que se tarda en que ocurra un determinado suceso. Por el nombre de la técnica parecería que se analizara el tiempo hasta la muerte (Análisis de supervivencia) pero, en realidad, puede analizarse cualquier otro suceso.

En Análisis de supervivencia la muestra consiste en el seguimiento de una serie de individuos desde el inicio del estudio hasta su final y, ante una situación de este tipo, es frecuente que se produzca la desaparición de alguno de esos individuos que entran en el estudio. También es posible que al entrar un individuo en el estudio, éste termine antes de que en ese individuo se produzca el suceso que se pretende detectar. Un ejemplo es el tiempo en que cominezan a aparecer determinados sintomas en pacientes sanos:

*Para modelar tablas de contingencia.Modelando tablas de contingencia (cruzar dos características y ver el número de casos que pertenecen a cada categoria)

qownnotes-media-uSdaBj

*Aproximaciones binomiales cuando n es muy grande y p muy pequeño, por ejemplo en epidemiología el tener o no tener una enfermedad sobre una población.

qownnotes-media-p14800

qownnotes-media-l14800

qownnotes-media-n19368

qownnotes-media-m19368

ESTUDIO ASINTÓTICO

Es el término para el análisis del comportamiento de las estadísticas cuando el tamaño de la muestra tiene un límite (normalmente infinito). Es la base de estudio frecuentista de probabilidades.

Un estimador es consistente si converge con lo que queremos estimar. La ley de números largos (LLN) dice que la media de la muestra, la varianza y la desviación estandar de iid muestras es consistente con la media de la población, la varianza y la desviación estándar respectivamente.

Teorema central del limite (CLT) : dice que la distribución de las medias siguen una distribución normal estándar centrada en el origen (ya que estamos restando la media de la población), con una desviación estándar igual al error estándar de la media, pero debe ser para un número de repeticiones suficiente como para que se cumpla.

\(X_1,\ldots,X_n\) es una colección de variables aleatorias iid con media \(\mu\) y varianza \(\sigma^2\) \(\bar X_n\) es la media de la muestra Entonces \(\frac{\bar X_n - \mu}{\sigma / \sqrt{n}}\) tiene una distribución normal para largos \(n\).

Por ejemplo, si tiramos un dado muchas veces y hacemos n observaciones

\(\frac{\bar X_n - 3.5}{1.71/\sqrt{n}}\)

Se ve claramente la dierencia de hacer 1,2 y 6 tiradas 10000 veces. La curva normal estándar lógicamente va a estar centrada en el origen ya que le estamos restando la media.

par(mfrow = c(1, 3))
for (n in c(1, 2, 6)){
  temp <- matrix(sample(1 : 6, n * 10000, replace = TRUE), ncol = n)
  temp <- apply(temp, 1, mean)
  temp <- (temp - 3.5) / (1.71 / sqrt(n))
  dty <- density(temp)
  plot(dty$x, dty$y, xlab = "", ylab = "density", type = "n", xlim = c(-3, 3), ylim = c(0, .5))
  title(paste("Media de la muestra con ", n, "observaciones"))
  lines(seq(-3, 3, length = 100), dnorm(seq(-3, 3, length = 100)), col = grey(.8), lwd = 3)
  lines(dty$x, dty$y, lwd = 2)
}

Al seguir una distribución normal, podemos usar los cuantiles para saber con qué probabilidad la media obtenida estará incluida en ese intervalo, para una distrióbucin normal sabemos \[\bar X_n \pm z_{1-\alpha/2}\sigma / \sqrt{n}\] y lo llamaremos intervalo de confianza.

En este ejemplo de las alturas de los hijos, podemos dar un intervalo de confianza en el cuál con un 95% de posibilidades la media obtenida en la muestra está contenida en la media de la población aplicando esta ecuación

library(UsingR);data(father.son);
x <- father.son$sheight
#Intervalo de confianza calculado a mano
(mean(x) + c(-1, 1) * qnorm(.975) * sd(x) / sqrt(length(x))) / 12

## [1] 5.709670 5.737674

Esta misma operación la podemos hacer con Bernuilli para proporciones de éxito o fallo:

\[ \hat p \pm z_{1 - \alpha/2} \sqrt{\frac{p(1 - p)}{n}} \] Sabiendo que la moneda es imparcial, p = 1/2, p(1-p)=1/4 y sqrt(p(1-p)) = 1/2

Aplicando le intervalor Wald del 95%

\(p(1-p) \leq 1/4\) for \(0 \leq p \leq 1\) \(\alpha = .05\) so that \(z_{1 -\alpha/2} = 1.96 \approx 2\) \[ 2 \sqrt{\frac{p(1 - p)}{n}} \leq 2 \sqrt{\frac{1}{4n}} = \frac{1}{\sqrt{n}} \] Entonces \(\hat p \pm \frac{1}{\sqrt{n}}\) es un buen intervalo de confianza rápido para \(p\)

Un ejemplo, es que en un sondeo de 100 personas, 56 votan sí, tengo la votación ganada? 1/sqrt(100)=.1 esto nos da un intervalo de confianza de Wald (95%) de entre (0.46, 0.66) así que no tenemos nada claro.

(0.56 + c(-1, 1) * qnorm(.975) * sqrt(0.56*0.44/100))

## [1] 0.4627099 0.6572901

binom.test(56,100)$conf.int

## [1] 0.4571875 0.6591640
## attr(,"conf.level")
## [1] 0.95

Si el proceso no fuera parcial (p=0.9 por ejemplo) el CTL nos garantiza únicamente que es una buena aproximación en el límite, pero se puede observar que si hacemos la misma prueba que con p = 0.5 tardan muchas más iteraciones en converger hacia una normal estándar.

Hay veces que el intervalo de wald no se comporta de manera correcta si el número de repeticiones no es el adecuado, y para evitarlo se puede aplicar la corrección de Agresti/Coull añadiendo dos aciertos y dos errores $p = X_n +2 / n+4 $

Por último si se trata de una distribución de Poisson por ejemplo para conteos:

Vemos que damos un un error diagnóstico 5 veces cada 94.32 días, cuál es el porcentaje de fallo diario con un intervalo de confianza dentro del 95%?

\(X \sim Poisson(\lambda t)\) Estimate \(\hat \lambda = X/t\) \(Var(\hat \lambda) = \lambda / t\)

\[ \frac{\hat \lambda - \lambda}{\sqrt{\hat \lambda / t}} = \frac{X - t \lambda}{\sqrt{X}} \rightarrow N(0,1) \]

x <- 5
t <- 94.32
lambda <- x/t
round(lambda + c(-1, 1) * qnorm(.975) * sqrt(lambda/t),3 )

## [1] 0.007 0.099

poisson.test(5,94.32)$conf.int

## [1] 0.01721254 0.12371005
## attr(,"conf.level")
## [1] 0.95

En este caso podemos observar que el intervalo calculado es más conservador ya que t no es lo suficientemente grande.

qownnotes-media-bG7876

Hay veces que el número de repeticiones no es suficiente en distribuciones discretas binomiales y de Poisson y se puede hacer una aproximación del valor estimado sumando +2 al numerador y +4 al denominador (añadiendo 2 exitos y 2 fallos)

qownnotes-media-jY7876

qownnotes-media-ZUuUAB

INTERVALOS DE CONFIANZA

Resumiendolo que hemos visto anteriormente, la media de los valores medios se aproxima a una normal estandar, con media \(\mu\) y sd \({\sigma / \sqrt{n}}\).

Por tanto, al seguir una normal estandar, el desplazamiento desde la media de 2 unidades del error estándar nos dará un intervalo en el que tendremos la probabilidad del 95% de contener el valor medio. A este intervalo lo llamamos ZQuantile X SE

\(\bar X_n \pm 2\sigma / \sqrt{n}\) es el llamado intervalo 95% de \(\mu\)

En la vida real, las muestras con las que trabajamos no suelen ser tan grandes, por lo que necesitamos una adaptación, utilizando la distribución T Student de Gosset. surge del problema de estimar la media de una población normalmente distribuida cuando el tamaño de la muestra es pequeño. Está también centrada en 0. Si no fuera simétrica en 0 (skewed distr) habría que calculalo de otra manera o coger otro estimador como la mediana

qownnotes-media-WYnisx

Aparece de manera natural al realizar la prueba t de Student para la determinación de las diferencias entre dos medias muestrales y para la construcción del intervalo de confianza para la diferencia entre las medias de dos poblaciones cuando se desconoce la desviación típica de una población y ésta debe ser estimada a partir de los datos de una muestra.

qownnotes-media-MIqXgZ

En adelante vamos a ver pequeñas muestras de datos y en este caso tendremos que usar la distribución T de Gosset y intervalos de confianza T, la única diferencia con la ecuación anterior es que usaremos el TQuantile

Cuando una distribución T creece y crece en muestras se acaba convirtiendo en distribución Z.

La distribución T tiene los extremos más pequños que la normal.A diferencia de la dist estandar, que tiene 2 parametros (media y varianza). La T de gosset tiene sólo un parámetro llamado grados de libertad t = n-1, cuanto más crezca este número más se parecerá a una standar normal, ya que si n es muy grande la diferencia entre una estandar normal, que es dividida por sigma, y la T de student que es dividida por la varianza, se hace insignificativa, es algo similar a lo que nos pasaba al estudiar la varianza de las muestras, que cuanto mayor era el número de las muestras más se acerca al teórico de la población. En este caso para calcular T con n-1 grados de libertad:

qownnotes-media-KN1524

El valor de esta t nos dará un intervalo de confianza, para el cuál dependiendo del número de individuos, en 95% de los casos ese valor está incluido, en el caso de que quisiéramos usar ese quantil qt(0.975,n)

Dependiendo de los grados de confianza, la diferencia de las distribuciones será mayor, en este caso los grados de confianza son muy amplios, lo que hace que sea muy similar

qownnotes-media-ASAszo

La distribución T es muy útil cuando queremos comparar cómo varía un grupo frente a otro, pero siempre teniendo en cuenta que deben ser datos iid normales, centrados en su media

Vamos a coger dos grupos y vamos a estudiar la diferencia de variabilidad entre ellos, primero estando relaccionados (que los individuos sean los mismos, y queremos ver su evolución) y luego siendo grupos independientes (Son grupos de personas totalmente distintas) .

Esa diferencia de variabilidad estará encerrada en un intervalo de confianza y se suele utilizar el intervalo T para compararlos

Los resultados del estudio por tanto varían dependiendo de si se emparejan los resultados porque están relaccionados o son independientes: Los datos son un grupo de 10 pacientes a los que se les ha tratado con 2 drogas distintas y se intenta comparar los resultados

data(sleep)
sleep

##    extra group ID
## 1    0.7     1  1
## 2   -1.6     1  2
## 3   -0.2     1  3
## 4   -1.2     1  4
## 5   -0.1     1  5
## 6    3.4     1  6
## 7    3.7     1  7
## 8    0.8     1  8
## 9    0.0     1  9
## 10   2.0     1 10
## 11   1.9     2  1
## 12   0.8     2  2
## 13   1.1     2  3
## 14   0.1     2  4
## 15  -0.1     2  5
## 16   4.4     2  6
## 17   5.5     2  7
## 18   1.6     2  8
## 19   4.6     2  9
## 20   3.4     2 10

En las siguientes gráficas se ven cada uno de los valores de ambos resultados, estando la leyenda la correspondencia entre el color de la raya que une a los individuos y el número de individuo que es en el grupo

DEPENDIENTES (paired=TRUE)

## [1] -3.363874  0.203874

qownnotes-media-JY1524

qownnotes-media-rV1524

Si queremos calcular el intevalo de confianza de los resultados emparejados, creamos un vector nuevo con la diferencia de los resultados de ambos grupos difference <- g2-g1

y se calcula la media y la desviación estandar mn <- mean(difference) s <- sd(difference) Los calculos se pueden hacer a mano mn + c(1,-1)qt(.975)(s/sqrt(10)) o con la función t.test()

qownnotes-media-qs1524

qownnotes-media-ay1524

Esto lo que nos dice, es que con 95% de posibilidades ( la funcion qt(.975) nos da la t para el 95% con 9 grados de libertad) la diferencia media de los efectos entre ambas drogas es entre 0.7 y 2.46 horas de sueño

data(sleep)
x1 <- sleep$extra[sleep$group == 1]
x2 <- sleep$extra[sleep$group == 2]
n1 <- length(x1)
n2 <- length(x2)
sp <- sqrt( ((n1 - 1) * sd(x1)^2 + (n2-1) * sd(x2)^2) / (n1 + n2-2))
md <- mean(x1) - mean(x2)
semd <- sp * sqrt(1 / n1 + 1/n2)
md + c(-1, 1) * qt(.975, n1 + n2 - 2) * semd

## [1] -3.363874  0.203874

t.test(x1, x2, paired = FALSE, var.equal = TRUE)$conf

## [1] -3.363874  0.203874
## attr(,"conf.level")
## [1] 0.95

t.test(x1, x2, paired = TRUE)$conf

## [1] -2.4598858 -0.7001142
## attr(,"conf.level")
## [1] 0.95

INDEPENDIENTES (PAIRED= FALSE)

En este caso tenemos dos grupos distintos de personas, unos que tomaron el medicamento y otro que tomaron placebo y vamos a medir su presión arterial. No se pueden emparejar porque son personas distintas y los grupos pueden tener incluso tamaño diferente.

qownnotes-media-Fz1524

El objetivo de nuevo es encontrar el intervalo de confianza de la diferencia media de la presión arterial media entre ambos grupos con un 95% de posibilidad. Para hacer el cálculo de los intervalos de confianza vamos a trabajar restando la media de ambos conjuntos, el quantile se calcula en base a ambos conjuntos, ya que hemos calculado su media por separado, así que hay que restar 2 a la suma de ambos para obtener el grado de confianza de todo.

qownnotes-media-aV1524

La varianza común es más dificil de calcular, haremos una estimación media entre ambas varianzas aplicando pesos en cada una de ellas dependiento del tamaño de cada grupo diviendolo por el tamaño total. Se ve más claro en este ejemplo extraído del libro Fundamentals of Biostatistics. Primero haremos los cálculos suponiendo que las varianzas son equivalentes. Se puede observar que el intervalo de confianza sobre difrencia media de ambas presiones arteriales medias contiene el valor 0, por lo que no se puede descartar que las medias de ambos grupos sean iguales.

qownnotes-media-Tn1524

Aquí se puede apreciar el error de cálculo si tratamos el ejemplo anterior como grupos independientes

qownnotes-media-kP1524

Se puede apreciar en el cálculo de ambos casos que cuando los tratamos de manera distinta un cálculo contiene el 0 y el otro no lo contiene, lo que significa que cuando los grupos están emparejados tienen muchas más medias distintas, de manera opuesta a cuando no los emparejamos y los tratamos como independientes en que se contiene el valor 0 y por tanto las diferencias en las medias son menos distintas ya que se tratan todos los datos como un solo grupo.

La función t.test tiene otro parámetro con el que se puede jugar que es si la varianza de los datos es igual en ambos conjuntos. En este ejemplo, se intenta hacer un intervalo de confianza entre la ganancia de peso de pollos con 4 dietas, comparándolas de 2 en 2, en este ejmplo la 1 y la 4, cuya diferencia de variablidad es muy evidente a nivel gráfico.

qownnotes-media-Bt1524

qownnotes-media-mA1524

vemos que los resultados de las estimaciones varían dependiendo del tratamiento de la relacción entre las varianzas de los grupos.

qownnotes-media-ll1524

DISTINTA VARIANZA ENTRE GRUPOS TEST DE WELCH

En el caso de que la varianza de los grupos no sea la misma o no estemos seguros, En caso de duda, se debe utilizar varianza no es igual, para calcularala

qownnotes-media-Mp1524

qownnotes-media-Rh1524

CHI SQUARE TEST (DATOS CATEGÓRICOS)

Este test se utiliza para saber si dos variables son independientes, es decir, si la aparición de una está relaccionada (ojo no es la causa, sólo si están relacconada) con la otra. Es apropiado utilizar este método en las siguientes situaciones:

• El muestreo es aleatorio
• las variables son categoricas
• si los datos recogidos en la tabla de contingencia tienen en las frecuencias esperadas al menos 5 casos.

Este método consiste en 4 pasos:

Establecer la hipótesis: Ho-> Variable A y variable B son independientes Ha-> Variable A y variable B no son independientes . La hipótesis alternativa por tanto es que sabiendo el valor A nos puede ayudar a predecir el valor de B Formular el plan de análisis Este paso describe cómo usar muestras de los datos para aceptar o rechazar la hipótesis nula: Nivel de significancia elegir un nivel apropiado entre 0 y 1 Test a utilizar en este caso chi-cuadrado.

Analizar la muestra de los datos Grados de libertad DF = (r-1)(c-1) donde r son el nº de niveles de A y c el nº de niveles de B Frecuencia esperada = Son computadas por separado para cada nivel de cada variable en cada nivel de la otra variable Er,c = (nr*nc)/n Donde Er,c es la frecuencia esperada para el nivel r de la variable A y nivel c de la B nr es el número total de observaciones del nivel r en A nc es el número total de observaciones del nivel c en B n es el número total de la muestra Test estadístico = x2 = sum[(Or,c - Er,c)2/Er,c] Donde Or,c es la frecuencia del nivel r de la variable A y nivel C de la B Er,c es la frecuencia esperada del nivel r de la variable A y el nivel c de la B P-value = es la probabilidad de observar una muestra de manera aleatoria al menos tan extrema como la nuestra. Para poder conocer esta probabilidad, se debe calcular si la distribución pertenece a una chi-cuadrado, si es así, utilizar los grados de libertad anteriores y el test estadístico para obtener el pvalue de la muestra de datos.

CONTRASTE DE HIPÓTESIS (Hypothesis testing) DATOS CON DISTRIBUCIONES NORMALES

Un contraste de hipótesis (también denominado test de hipótesis o prueba de significación) es un procedimiento para juzgar si una propiedad que se supone en una población estadística es compatible con lo observado en una muestra de dicha población.La aplicación de cálculos probabilísticos permite determinar a partir de qué valor debemos rechazar la hipótesis garantizando que la probabilidad de cometer un error es un valor conocido a priori.

En el contraste de hipótesis tenemos sempre dos resultados:

Una hipotesis nula: H0 (de no efecto) Una hipotesis alternativa: H1 (de diferencia o efecto)

A partir de una muestra de la población en estudio, se extrae un estadístico (media muestral, varianza muestral,etc.) cuya distribución de probabilidad esté relacionada con la hipótesis en estudio y sea conocida. Se toma entonces como región de rechazo al conjunto de valores que es más improbable bajo la hipótesis, esto es, el conjunto de valores para el que rechazaremos la hipótesis nula si el valor del estadístico observado entra dentro de él.

Para probar una hipótesis vamos a comparar siempre con la hipótesis nula Ho, e invalidarla, ya que partimos de que la hipótesis nula es siempre cierta.

A continuación, vamos a ver cuatro ejemplos de cómo descartar H_0, con una distribución normal grande, con una distribución normal pequeña ambos con varianza equvalente, uno con varianza no equivalente y un ejemplo con distribición binomial.

Z-TEST

Vamos a trabajar con un conjunto de 100 pacientes con sobrepeso que tienen desorden del sueño y una media de 32 eventos/hora durante la noche con una desviación estándar de 10 eventos/hora. Se considera que desorden del sueño es 30 eventos/hora. Vamos a utilizar el Test-Z como test estadístico, para poder usarlo necesitamos la media , la varianza y la distribución.

Queremos invalidar H0 que dice que mu = 30 y validar H_a que dice que mu > 30 (32 en este caso concreto)

ERROR TIPO I –> EL investigador SI encuentra diferencia entre Ho y Ha cuando no la hay (Falso positivo) ALFA probabilidad de cometer ese error

ERROR TIPO II –> El investigador NO encuentra diferencia entre Ho y Ha cuando sí la hay (falso negativo) BETA probabilidad de cometer ese error

Error tipo I –> Rechazar Ho cuando es correcta esto es alfa (probabilidad de cometer este error es típicamente 5% llamada alfa) Tampoco queremos un alfa muy bajo ya que nunca rechazaríamos H_o aunque fuera falsa. Error tipo II –> Aceptar H_o cuando es falso Error tipo III –> Rechazar H_a cuando es correcta Error tipo IV –> Aceptar H_a cuando es falsa

qownnotes-media-amRzlL

El error estándar de la media de un conjunto es desviacion_estandar/sqrt(población), en este ejemplo 10/sqrt(100) = 1

H_o por tanto tiene un media normalmente distribuida,mu = 30 y varianza = 1 (estimada como el error estandar que es 1 al cuadrado)

Estamos buscando una constante C que hiciera que la probabilidad de que la media sea mayor que C es %5 Esto es P(X>C|H_o) = 5%

!(media/32517.png)

El área sin colorear se puede calcular como qnorm(.95)

Recordando que en una distribución normal el percentil 95 está a 1.645 desviaciones estándar de la media y teniendo una distribución normal de media 30 y 1 de desviación estandar N(30,1)

Rechazaríamos la hipotesis nula si la media fuera mayor que un valor (C) cuya probabilidad sea 5% C = 30 + 1x1.645 = 31.645

!(media/26308.png)

Rechazar la hipótesis nula cuando la media es >= 31.645 tiene una probabilidad de fallo del 5%, nuestra media es 32 que es mayor y cae dentro de ese área así que rechazamos H_o En este gráfico vemos las dos distribuciones, y cómo la media de H_a está dentro del área del 5% de H_o, lo que nos sirve para rechazar que nuestro resultado sea por azar en un 95% de los casos al menos.

!(media/26001.png)

El Z-score es lo mismo pero calculando cuantos errores estándar la media de la hipótesis es diferente de la media, esto se calcula computando la distancia entre las medias (32-30) y dividiendo por el error estandar de la media (s/sqrt(n)) que en este caso es 10/sqrt(100)

Z-score = 2, significa que está a 2 errores estandar de la media.

Un alumno A saca una puntuación de 85 en un examen cuyas puntuaciones tienen una media de 79 con una desviación típica de 8. Un alumno B saca 74 en un examen cuyas puntuaciones tienen una media de 70 y desviación típica de 5. ¿Cuál de los dos alumnos obtuvo una puntuación mejor? La respuesta, desde el punto de vista de la unidad tipificada, se obtiene así:

Las puntuaciones tipificadas de los alumnos A y B son respectivamente:

qownnotes-media-blJJQy

Así el alumno B lo hizo mejor que el A, aunque su puntuación de 74 es inferior a 85.

qownnotes-media-rn7296

La regla por tanto para descartar H_0 es (X-mu)*sqrt(n)/s > Z

T-TEST

Digamos que ahora bajamos de 100 ejemplos a 16, se convierte entonces en una distribución T con 15 grados de libertad. Aquí ya no podemos utlizar Zscore como test estadístico porque la distribución es normal, pero con pocos individuos (distribución T-student)

qownnotes-media-Hy7296

En este caso fallamos al rechazar Ho porque 0.8< 1.7531

En lugar de concentrar todo el error en uno de los extremos, se suele distribuir a ambos lados

qownnotes-media-sv7296

Normalmente todos estos cálculos no los tenemos que hacer a mano, ya que la función t.test() nos da toda esta información

qownnotes-media-nS7296

En este ejemplo, queremos saber si la media de altura de los hijos es equivalente a la media de altura de los padres. Cogemos los datos de las alturas de padres e hijos, y restamos la altura de un padre a la de su hijo, y vamos a ver si la diferencia (haciendo la resta se asume que están emparejados) es 0 o no es 0 con t.test (tambien se podrían pasar ambos arrays a la función y como tercer parametro paired.equals=TRUE)

Ambos t(cuantas desviaciones estandar está la media del conjunto del de H_o) y df (grados de libertad) son muy grandes con lo que no hace falta calcularlo.

Otro ejemplo

qownnotes-media-cy7296

library(effsize)
datos <- data.frame(corredor = c(1:10),
                    antes = c(12.9, 13.5, 12.8, 15.6, 17.2, 19.2, 12.6, 15.3,
                              14.4, 11.3),
                    despues = c(12.7, 13.6, 12.0, 15.2, 16.8, 20.0, 12.0, 15.9,
                                16.0, 11.1))
cohen.d(d = datos$antes, f = datos$despues, paired = TRUE)

##
## Cohen's d
##
## d estimate: -0.06745406 (negligible)
## 95 percent confidence interval:
##     lower     upper
## -1.007282  0.872374

CONTRASTE DE HIPÓTESIS (Hypothesis testing) NO NORMALES

A la hora de hacer contraste de hipótesis, podemos hablar de un contraste cuantitativo (cuán diferentes es una ditribución de algo medido en datos CONTINUOS, por ejemplo días de estancia, respecto de otra, ya sea su distribución supuesta u otra) o un contraste cualitativo (basado en el conteo de ocurrencias de variables categóricas DISCRETOS, por ejemplo tipos de altas al mes, si dos conteos tienen diferencias significativas, ya sea frente a su distribución supuesta u otra).

Otro factor fundamental a la hora de hacer el contraste será si conocemos su distribución supuesta (Test paramétricos, porque usaremos sus parámetros poblacionales media y varianza) o si no podemos definirla (No paramétricos).

Para hacer esta primera valoración, miraremos si cumple los criterios de normalidad y si sus varianzas son equivalentes.

Test de Normalidad Shapiro-Wilk Este test se emplea para contrastar normalidad cuando el tamaño de la muestra es menor de 50. Para muestras grandes es equivalente al test de kolmogorov-Smirnov. Uno de los supuestos que debemos de tratar en este tipo de prueba es el de distribucion normal. Utilizaremos para ello el estadistico de Shapiro-Wilk, existen otros estadisticos que pueden tambien utilizarse para probar la distribucion normal.El problema que nos encontramos es que para más de 5000 elementos este método ineficaz. Test de Kolmogorov-Smirnov Este es otro tipo de prueba que se utiliza para comparar la normalidad de un conjunto de datos. La prueba determina si los datos provienen de una misma distribucion. Puede ser utilizado tambien para comparar si un conjunto de datos se ajustan a una distribucion (en nuestro caso normal).

Cuando las muestras son pequenas la prueba de Kolmogorov es una prueba menos poderosa para detectar diferencias significativas.

El test de Kolmogorov-Smirnov permite estudiar si una muestra procede de una población con una determinada distribución (media y desviación típica), no está limitado únicamente a la distribución normal. Se ejecuta con la función ks.test().

ks.test(x = mtcars$mpg,"pnorm", mean(mtcars$mpg), sd(mtcars$mpg))

## Warning in ks.test(x = mtcars$mpg, "pnorm", mean(mtcars$mpg),
## sd(mtcars$mpg)): ties should not be present for the Kolmogorov-Smirnov test

##
##  One-sample Kolmogorov-Smirnov test
##
## data:  mtcars$mpg
## D = 0.1263, p-value = 0.687
## alternative hypothesis: two-sided

A pesar de que continuamente se alude al test Kolmogorov-Smirnov como un test válido para contrastar la normalidad, esto no es del todo cierto. El Kolmogorov-Smirnov asume que se conoce la media y varianza poblacional, lo que en la mayoría de los casos no es posible. Esto hace que el test sea muy conservador y poco potente. Para solventar este problema, se desarrolló una modificación del Kolmogorov-Smirnov conocida como test Lilliefors. El test Lilliefors asume que la media y varianza son desconocidas, estando especialmente desarrollado para contrastar la normalidad. Es la alternativa al test de Shapiro-Wilk cuando el número de observaciones es mayor de 50. La función lillie.test() del paquete nortest permite aplicarlo.

library("nortest")
lillie.test(x = mtcars$mpg)

##
##  Lilliefors (Kolmogorov-Smirnov) normality test
##
## data:  mtcars$mpg
## D = 0.1263, p-value = 0.2171

Respecto a la prueba de normalidad de Shapiro, esta resulta ser mas sensitiva que Kolmogorov. Es probable que la primera muestre que no existe normalidad y la segunda que si.

P1<-c(111,128,111,117,119,90,115,118,109,115,118,113,92,79,75,93,86,99,118)
P2<-c(88,118,100,107,119,87,105,121,111,117,108,103,96,82,72,90,91,92,115)
shapiro.test(P1)

##
##  Shapiro-Wilk normality test
##
## data:  P1
## W = 0.89668, p-value = 0.04239

mean<-mean(P1)
sd<-sd(P1)
ks.test(P1, "pnorm", mean=mean(P1), sd=sd(P1))

## Warning in ks.test(P1, "pnorm", mean = mean(P1), sd = sd(P1)): ties should
## not be present for the Kolmogorov-Smirnov test

##
##  One-sample Kolmogorov-Smirnov test
##
## data:  P1
## D = 0.21985, p-value = 0.3174
## alternative hypothesis: two-sided

Como observa, en este caso Shapiro muestra que no existe normalidad de datos (P1), mientras que Kolmogorov si lo muestra, en el mismo vector, demostrando asi la sensibilidad de cada una de las pruebas.

ks.test(P1,P2)

## Warning in ks.test(P1, P2): cannot compute exact p-value with ties

##
##  Two-sample Kolmogorov-Smirnov test
##
## data:  P1 and P2
## D = 0.31579, p-value = 0.2997
## alternative hypothesis: two-sided

Test de Jarque-Bera El test de Jarque-Bera no requiere estimaciones de los parámetros que caracterizan la normal. El estadístico de Jarque-Bera cuantifica que tanto se desvían los coeficientes de asimetría y curtosis de los esperados en una distribución normal.

Los datos presentan normalidad. Recuerde siempre visualizar de forma grafica los datos, tanto por su función de densidad (PDF) como por su función de distribución (CDF)

par(mfrow=c(2,2))
hist(P1)
hist(P2)
plot(ecdf(P1),xlim=range(c(P1,P2)),main="Distribucion Empirica Acumulada\n(P1,P2)")
plot(ecdf(P2), add=TRUE, lty="dashed", col="red")

par(mfrow=c(2,2))
qqnorm(P1,col="red")
qqline(P1)
qqnorm(P2,col="red")
qqline(P2)
qqplot(P1,P2, main="Q-Q Plot de P1 y P2")

Si las gráficas tienen la misma distribución, el gárfico Q-Q entre grupos mostrará una diagonal perfecta, cuanto más se ciñan los puntos a esa diagonal perfecta más similar será distribución. El diagrama Q-Q consiste en comparar los cuantiles de la distribución observada con los cuantiles teóricos de una distribución normal con la misma media y desviación estándar que los datos. Cuanto más se aproximen los datos a una normal, más alineados están los puntos entorno a la recta.

Para la representación gráfica se puede hacer también uso de esta librería:

library(PASWR)

## Loading required package: e1071

##
## Attaching package: 'e1071'

## The following object is masked from 'package:Hmisc':
##
##     impute

EDA(P1)

## [1] "P1"

## Size (n)  Missing  Minimum   1st Qu     Mean   Median   TrMean   3rd Qu
##   19.000    0.000   75.000   92.000  105.579  111.000  105.579  118.000
##     Max.   Stdev.     Var.  SE Mean   I.Q.R.    Range Kurtosis Skewness
##  128.000   15.334  235.146    3.518   26.000   53.000   -1.064   -0.585
## SW p-val
##    0.042

Test de la varianza F de Fisher El F-test, también conocido como contraste de la razón de varianzas, contrasta la hipótesis nula de que dos poblaciones normales tienen la misma varianza. Es muy potente, detecta diferencias muy sutiles, pero es muy sensible a violaciones de la normalidad de las poblaciones

Se utiliza para saber si las varianzas de dos distribuciones son iguales. Como siempre: H0 las varianzas son iguales H1 las varianzas varían significativamente

var.test(P1,P2)

##
##  F test to compare two variances
##
## data:  P1 and P2
## F = 1.1874, num df = 18, denom df = 18, p-value = 0.7195
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.4574801 3.0820885
## sample estimates:
## ratio of variances
##           1.187432

Test de Bartlet Permite contrastar la igualdad de varianza en 2 o más poblaciones sin necesidad de que el tamaño de los grupos sea el mismo. Es más sensible que el test de Levene a la falta de normalidad, pero si se está seguro de que los datos provienen de una distribución normal, es la mejor opción.

data("iris")
a <- iris[iris$Species == "versicolor", "Petal.Length"]
b <- iris[iris$Species == "virginica", "Petal.Length"]
bartlett.test(list(a,b))

##
##  Bartlett test of homogeneity of variances
##
## data:  list(a, b)
## Bartlett's K-squared = 1.249, df = 1, p-value = 0.2637

Test de varianza Levene El test de Levene se caracteriza, además de por poder comparar 2 o más poblaciones, por permitir elegir entre diferentes estadísticos de centralidad :mediana (por defecto), media, media truncada. Esto es importante a la hora de contrastar la homocedasticidad dependiendo de si los grupos se distribuyen de forma normal o no.

library(car)

## Loading required package: carData

##
## Attaching package: 'carData'

## The following object is masked from 'package:PASWR':
##
##     Wool

data("iris")
iris <- dplyr::filter(.data = iris, Species %in% c("versicolor", "virginica"))
leveneTest(y = iris$Petal.Length, group = iris$Species, center = "median")

## Levene's Test for Homogeneity of Variance (center = "median")
##       Df F value Pr(>F)
## group  1  1.0674 0.3041
##       98

Este test se debe aplicar cuando se quieran comparar más de dos grupos como en un estudio ANOVA. Un p-value por debajo de 0.05 nos dice que las varianzas de los datos nos iguales, por tanto no se pueden usar test paramétricos como ANOVA.

library(car)
size<-c(25,22,28,24,26,24,22,21,23,25,26,30,25,24,21,27,28,23,25,24,20,22,24,23,22,24,20,19,21,22)
location<-c(rep("ForestA",10), rep("ForestB",10), rep("ForestC",10))
my.dataframe<-data.frame(size,location)
car::leveneTest(size~location, data=my.dataframe)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.4919 0.6168
##       27

Este test nos da un p-value muy alto lo que significa que no hay diferencias estadísticamente significativas. Se puede el parámetro center = mean para que use la media como centro en lugar de la mediana por defecto, que hace el test más robusto y menos sensible a outliers. En este caso el test se llama Brown-Forsythe

Test de varianza Fligner-Killeen

Si los datos no superasen el test de normalidad, deberíamos aplicar este otro método ya que es una prueba menos sensible a que los datos sean normales o cuando tengamos outliers no resueltos

#Debemos hacer también el test de Fligner ya que es una prueba menos sensible a que los datos sean normales
fligner.test(x = list(P1,P2))

##
##  Fligner-Killeen test of homogeneity of variances
##
## data:  list(P1, P2)
## Fligner-Killeen:med chi-squared = 0.04751, df = 1, p-value =
## 0.8275

#H0= Los datos tienen varianzas iguales H1: Los datos presentan varianzas distintas.

Datos emparejados

Lo último que podemos necesitar saber es si los sujetos de los grupos están o no emparejados. Si esto no lo sabemos a priori podemos hacer un análisis de correlación, si los datos tinen un alto nivel de correlación es un indicio de que estén emparejados.

cor(P1, P2)

## [1] 0.8838955

#H0= Los datos tienen varianzas iguales H1: Los datos presentan varianzas distintas.

Podemos comprobar gráficamente que están muy relacionadas

plot(P1, P2, pch=20, main = "P1 y P2")

a = 24 + 10*rbeta(150, 1.1, 1.1)  # generate fake data
b = 24 + 10*rbeta(150, 1.1, 1.1)
err = rnorm(150, 0, .1);  aa = a + err
plot(a, b, pch=20, main="Independent A and B");  plot(a, aa, pch=20, main="Paired A and AA")

Commo hemos podido comprobar, a y b están relaccionados a pesar de tener la misma distribución, mientras que P1 y P2 no parecen tener igual distribución aunque sí hemos visto que estaban relaccionadas:

qqplot(a,b,main="Q-QPlot A y B");qqplot(P1,P2,main="Q-QPlot P1 y P1")

Con toda esta infromación que hemos recabado, ya podemos confirmar si los datos son o no normales, si tienen o no distinta varianza y si se deben o no emparejar los resultados de los grupos o no.

muestraX <- c( 1.1, 3.4, 4.3, 2.1, 7.0 , 2.5 )
muestraY <- c( 7.0, 8.0, 3.0, 5.0, 6.2 , 4.4 )
fligner.test(x = list(muestraX,muestraY))

##
##  Fligner-Killeen test of homogeneity of variances
##
## data:  list(muestraX, muestraY)
## Fligner-Killeen:med chi-squared = 0.07201, df = 1, p-value =
## 0.7884

#No hay evidencias en contra de la igualdad de varianzas.
wilcox.test(x = muestraX, y = muestraY, alternative = "two.sided", mu = 0,
            paired = FALSE, conf.int = 0.95)

## Warning in wilcox.test.default(x = muestraX, y = muestraY, alternative =
## "two.sided", : cannot compute exact p-value with ties

## Warning in wilcox.test.default(x = muestraX, y = muestraY, alternative =
## "two.sided", : cannot compute exact confidence intervals with ties

##
##  Wilcoxon rank sum test with continuity correction
##
## data:  muestraX and muestraY
## W = 6.5, p-value = 0.07765
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -4.9000274  0.4000139
## sample estimates:
## difference in location
##              -2.390035

require(coin)

## Loading required package: coin

# La función wilcox.test() del paquete coin requiere pasarle los argumentos en
# forma de función (~), por lo que los datos tienen que estar almacenados en
# forma de data frame.
datos <- data.frame(grupo = rep(c("A", "B"), c(6, 6)),
                    valores = c(muestraX, muestraY))
wilcox_test(valores ~ grupo, data = datos, distribution = "exact", conf.int=0.95)

##
##  Exact Wilcoxon-Mann-Whitney Test
##
## data:  valores by grupo (A, B)
## Z = -1.8447, p-value = 0.06926
## alternative hypothesis: true mu is not equal to 0
## 95 percent confidence interval:
##  -4.9  0.4
## sample estimates:
## difference in location
##                   -2.4

n1 <- length(muestraX)
n2 <- length(muestraY)
# La función devuelve el valor de Z, por lo que se puede calcular fácilmente el
# tamaño del efecto
tamanyo_efecto <- 1.845/sqrt(n1 + n2)
tamanyo_efecto

## [1] 0.5326056

antes       <- c( 2, 5, 4, 6, 1, 3 )
despues     <- c( 5, 6, 2, 7, 1, 6 )

wilcox.test(x = antes, y = despues, alternative = "two.sided", mu = 0, paired = TRUE)

## Warning in wilcox.test.default(x = antes, y = despues, alternative =
## "two.sided", : cannot compute exact p-value with ties

## Warning in wilcox.test.default(x = antes, y = despues, alternative =
## "two.sided", : cannot compute exact p-value with zeroes

##
##  Wilcoxon signed rank test with continuity correction
##
## data:  antes and despues
## V = 3, p-value = 0.2763
## alternative hypothesis: true location shift is not equal to 0

require(coin)
# La función wilcoxsigned_test() del paquete coin requiere pasarle los
# argumentos en forma de función (~), por lo que los datos tienen que estar
# almacenados en forma de data frame.

datos <- data.frame(antes = antes, despues = despues)
wilcoxsign_test(antes ~ despues, data = datos, distribution = "exact")

##
##  Exact Wilcoxon-Pratt Signed-Rank Test
##
## data:  y by x (pos, neg)
##   stratified by block
## Z = -1.272, p-value = 0.25
## alternative hypothesis: true mu is not equal to 0

x <- seq(0, 5, length = 1000)
y <- dlnorm(x = x, meanlog =  0.5, sdlog = 1.7)
plot(x, y, type = "l", lty = 1, xlab = "x", col = "blue", ylab = "Densidad",
     main = "Distribuciones lognormal y normal con misma mediana", xlim=c(-5, 5))
x_2 <- seq(-5, 5, length = 1000)
y_2 <- dnorm(x = x_2, mean = 1.648721, sd = 1)
lines(x_2, y_2, col = "red")
legend("topright",
       legend = c("mediana lognormal = 0.5", "mediana normal = 0.5"),
       col = c("blue", "red"), lty = 1,cex = 0.7)

set.seed(888)
muestra_a <- rlnorm(n = 20, meanlog =  0.5, sdlog = 1.7)
muestra_b <- rnorm(n = 20, mean = 1.648721, sd = 1)
par(mfrow = c(1,2))
hist(muestra_a, breaks = 10, main = "", col = "blue")
hist(muestra_b, breaks = 10, main = "", col = "firebrick")

wilcox.test(muestra_a, muestra_b, paired = FALSE)

##
##  Wilcoxon rank sum test
##
## data:  muestra_a and muestra_b
## W = 196, p-value = 0.9254
## alternative hypothesis: true location shift is not equal to 0

set.seed(888)
muestra_a <- rlnorm(n = 1000, meanlog = 0.5, sdlog = 1.7)
muestra_b <- rnorm(n = 1000, mean = 1.648721, sd = 1)
wilcox.test(muestra_a, muestra_b, paired = FALSE)

##
##  Wilcoxon rank sum test with continuity correction
##
## data:  muestra_a and muestra_b
## W = 551900, p-value = 5.842e-05
## alternative hypothesis: true location shift is not equal to 0

CONTRASTE DE HIPÓTESIS DATOS CATEGÓRICOS

Los contrastes de hipótesis para variables cualitativas se realizan mediante test de frecuencia o proporciones. Dentro de esta categoría existen distintos tipos de test, la utilización de uno u otro depende de qué tipo de información se quiera obtener:

Test de distribución esperada o goodness of fit: Se emplean para comparar la distribución observada frente a una distribución esperada o teórica. Test de diferencia de frecuencias o test de independencia: Se emplean para estudiar si la frecuencia de observaciones es significativamente distinta entre dos o más grupos. En los test de goodness of fit solo hay una variable asociada a cada observación, mientras que en los test de independencia hay dos variables asociadas a cada observación. También se emplean distintos test dependiendo del tipo de datos (independientes o dependientes) con los que se vaya a trabajar. Las siguientes tablas muestran algunos de los más empleados.

Los test exactos calculan la probabilidad de obtener los resultados observados de forma directa generando todos los posibles escenarios y calculando la proporción en los que se cumple la condición estudiada (son test de permutaciones). Los test aproximados calculan primero un estadístico y luego emplean la distribución teórica de dicho estadístico para obtener la probabilidad de que adquiera valores iguales o más extremos.

Existe bastante controversia en cuanto a si se deben de utilizar test exactos o aproximados. En la era pre-computacional, los test exactos se complicaban mucho cuando el tamaño total de muestras aumentaba, sin embargo, por medio de la computación esta barrera se ha eliminado. Los test exactos son más precisos cuando el tamaño total de observaciones es bajo o alguno de los grupos tiene pocas observaciones, una vez alcanzado un número alto de observaciones las diferencias son mínimas. En el libro Handbook of Biological Statistics John H. McDonald se recomienda utilizar test exactos cuando el número total de observaciones es menor a 1000 o cuando, aunque el número total sea mayor a 1000, haya algún grupo cuyo número de eventos esperados sea pequeño (normalmente menor que 5). En el caso de aplicar test aproximados sobre tamaños pequeños se suelen emplear correcciones, las más frecuentes son la corrección de continuidad de Yate o la corrección de William.

Se puede considerar a los test basados el estadístico ??2 como una generalización del contraste de proporciones basado en la aproximación a la normal (Z-test de una proporción y Z-test de dos proporciones) cuando hay 2 o más variables cualitativas o alguna de ellas tiene 2 o más niveles. En aquellos casos en los que ambos test se pueden aplicar, el resultado de un Z-test y un test ??2 es equivalente. Esto es debido a que en la distribución chi-cuadrado con 1 grado de libertad el estadístico ??2 es igual al estadístico Z de una distribución normal, elevado al cuadrado. ##1 DISTRIBUCIÓN ESPERADA## ###TEST BINOMIAL EXACTO### El test binomial exacto se emplea para estudiar si la proporción de eventos verdaderos de una variable de tipo binomial se diferencia significativamente de la frecuencia teórica con la que se esperaría que apareciesen. En R puede realizarse un test binomial exacto mediante la función binom.test(), obteniendo tanto el p-value como el intervalo de confianza del valor de p en la población.

Un estudio sugiere que en una ciudad están naciendo más hombres que mujeres. Para determinar si es cierto se selecciona una muestra aleatoria de los niños nacidos durante los últimos 2 años y se identifica para cada uno el género. ¿Existen diferencias significativas para un nivel de significancia del 5%?

datos <- c("masculino", "masculino", "masculino", "masculino", "femenino",
           "masculino", "masculino", "masculino", "femenino", "femenino",
           "femenino", "femenino", "femenino", "femenino", "masculino",
           "masculino", "femenino", "masculino", "femenino", "masculino",
           "femenino", "masculino", "masculino", "masculino", "femenino",
           "masculino", "masculino", "femenino", "masculino", "femenino")
tabla <- table(datos)
tabla

## datos
##  femenino masculino
##        13        17

binom.test(x = tabla, alternative = "two.sided", conf.level = 0.95)

##
##  Exact binomial test
##
## data:  tabla
## number of successes = 13, number of trials = 30, p-value = 0.5847
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.2546075 0.6257265
## sample estimates:
## probability of success
##              0.4333333

El test binomial no encuentra diferencias significativas entre las observaciones y lo que cabría esperar si la verdadera probabilidad de nacimiento es del 50% para ambos sexos.

Más de 2 niveles

En los estudios genéticos de Mendel se establece que la descendencia de un determinado cruce de plantas da lugar a cuatro tipos de semillas (amarillas-lisas, verdes-lisas, amarillas-rugosas y verdes-rugosas) con una proporción teórica de 9:3:3:1. Para demostrar si es cierto se recoge una muestra de 556 semillas y se observa que las cantidades para cada grupo son de 315, 108, 101, 32 respectivamente. ¿Se corresponden las observaciones con el modelo esperado?

observaciones <- c(315, 108, 101, 32)
freq_teoricas <- c(9,3,3,1)

require(XNomial)

## Loading required package: XNomial

XNomial::xmulti(obs = observaciones, expr = freq_teoricas, statName = "Prob") 

##
## P value (Prob) = 0.9382

La función xmulti() también calcula el p-value a partir de modelos aproximados (explicados más adelante) y crea una representación gráfica que permite ver como de bien se ajusta la distribución chi-cuadrado a las probabilidades exactas.

xmulti(obs = observaciones, expr = freq_teoricas, detail = 3, histobins = TRUE)

##
## P value  (LLR)  =  0.9261
## P value (Prob)  =  0.9382
## P value (Chisq) =  0.9272
##
## Observed:  315 108 101 32
## Expected ratio:  9 3 3 1
## Total number of tables:  28956759

El p-value obtenido indica que las observaciones se asemejan mucho a las esperadas según la distribución de frecuencias del modelo teórico (Ho). No hay evidencias para rechazar la hipótesis nula.

CHI2 GOODNESS OF FIT

El test chi2 goodness of fit es la alternativa aproximada al test binomial (dos niveles) o multinomial (>2 niveles) normalmente utilizada cuando el número de observaciones o de grupos es demasiado alto para poder emplear los test exactos. Permite comparar la distribución de las observaciones en los diferentes niveles de una variable cualitativa con lo esperado acorde a una distribución hipotética Ho. Se trata por lo tanto de una expansión del Z-test para una proporción cuando la variable estudiada tiene dos o más niveles. Cuando la variable cualitativa tiene dos niveles, el test chi2 goodness of fit y el Z-test para una proporción son equivalentes.

El funcionamiento del test consiste en cuantificar y sumar las diferencias entre el número de eventos observadas en cada nivel con respecto al número esperado acorde a H0. El valor de este sumatorio se asigna a un estadístico llamado chi2 y se emplea su relación matemática con la distribución chi2 para estimar la probabilidad de obtener un valor igual o más extremo.

H0: Las observaciones siguen una distribución acorde a la distribución del modelo de referencia. Las desviaciones son debidas al azar.

Ha: Las observaciones no siguen una distribución acorde a la distribución del modelo de referencia. Las desviaciones son mayores a las esperadas por azar.

Un estudio intenta comparar si existe discriminación racial en la formación de jurados populares. Para ello se compara la distribución racial en el país con la distribución de los miembros de los últimos 500 jurados.

Ejemplo

chisq.test(x = c(1920, 347, 19, 84, 130), p = c(80.29, 12.06, 0.79, 2.92, 3.94) / 100)

##
##  Chi-squared test for given probabilities
##
## data:  c(1920, 347, 19, 84, 130)
## X-squared = 22.419, df = 4, p-value = 0.0001654

En caso de que no se cumpla el número mínimo de eventos por nivel, se recurre a la simulación:

chisq.test(x = c(1920, 347, 19, 84, 130),
           p = c(80.29, 12.06, 0.79, 2.92, 3.94)/100,
           simulate.p.value = TRUE, B = 5000)

##
##  Chi-squared test for given probabilities with simulated p-value
##  (based on 5000 replicates)
##
## data:  c(1920, 347, 19, 84, 130)
## X-squared = 22.419, df = NA, p-value = 0.0005999

Además del valor del p-value se recomienda incluir una representación gráfica de los valores observados y esperados.

require(ggplot2)
datos <- data.frame(raza=rep(c("blanca","negra","latina","asiática","otra"), 2),
                    tipo_valor = rep(c("observado", "esperado"), each = 5),
                    eventos = c(1920, 347, 19, 84, 130, 2007.25, 301.50,
                                19.75, 73.00, 98.50))
ggplot(data = datos, aes(x = raza, y = eventos)) +
geom_bar(stat = "identity", aes(fill = tipo_valor), position = "dodge") +
theme_bw()

Se confirma tanto por un pvalue muy bajo como gráficamente que no se está siguiendo la proporción esperada

2 COMPARACIÓN DE GRUPOS

TEST EXACTO DE FISHER

La prueba de Fisher es el test exacto utilizado cuando se quiere estudiar si existe asociación entre dos variables cualitativas, es decir, si las proporciones de una variable son diferentes dependiendo del valor que adquiera la otra variable. En la gran mayoría de casos, el test de Fisher se aplica para comparar dos variables categóricas con dos niveles cada una (tabla 2x2). Es posible utilizarlo con tablas 2xK niveles pero los requerimientos de cálculo son altos.

Ho: Las variables son independientes por lo que una variable no varía entre los distintos niveles de la otra variable.

Ha: Las variables son dependientes, una variable varía entre los distintos niveles de la otra variable.

El test de Fisher es más preciso que sus equivalentes aproximados (test chi-square de independencia o G???test de independencia) cuando el número de eventos esperado por nivel es pequeño. Se recomienda utilizarlo siempre que sea posible (tiempo de computación) aunque para observaciones totales >1000 los resultados de los test aproximados son muy parecidos.

Es importante tener en cuenta que el test de Fisher está diseñado para situaciones en las que las frecuencias marginales de filas y columnas (los totales de cada fila y columna) son fijas, se conocen de antemano. Esta condición es relevante en los experimentos biológicos ya que no es común poder cumplirla. Si esta condición no se satisface el test de Fisher deja de ser exacto, por lo general pasando a ser más conservativo. En varios artículos se menciona que el test de Barnard es más potente que el de Fisher cuando las frecuencias marginales no son fijas. También parece ser que aunque el test deja de ser exacto no significa que no se pueda aplicar.

Ejemplo de experimentos con y sin frecuencias marginales fijas:

Frecuencias marginales fijas: Supóngase que se quiere saber si la preferencia que tienen dos especies de pájaros (estorninos y gorriones) para refugiarse en casetas artificiales es diferente dependiendo del material de fabricación (madera o metal). Para ellos se disponen en una pajarera 5 casetas de metal y 5 de madera y se sueltan en el interior de la jaula 4 gorriones y 6 estorninos. En este experimento se sabe que las frecuencias marginales van a ser 5, 5, 4, 6 lo que no se sabe es como se van a distribuir las observaciones dentro de la tabla.

Tabla de métodos

Dado que el test de Fisher contrasta si las variables están relacionadas, al tamaño del efecto se le conoce como fuerza de asociación. Existen múltiples medidas de asociación, entre las que destacan phi o Cramer???s V. Los límites empleados para su clasificación son:

pequeño: 0.1 mediano: 0.3 grande: 0.5 En R se pueden calcular mediante la función assocstats() del paquete vcd.

Se quiere estudiar si la reacción alérgica a un compuesto y una determinada mutación en un gen están relacionados. Para ello se realiza un test alérgico sobre un grupo de individuos seleccionados al azar y se genotipa el estado del gen de interés ¿Existe un diferencia significativa en la incidencia de la mutación entre los alérgicos y no alérgicos?

datos <- data.frame( sujeto = c("No alérgico", "No alérgico", "No alérgico",
                                "No alérgico", "alérgico", "No alérgico",
                                "No alérgico", "alérgico", "alérgico",
                                "No alérgico", "alérgico", "alérgico",
                                "alérgico", "alérgico", "alérgico",
                                "No alérgico", "No alérgico", "No alérgico",
                                "No alérgico", "alérgico", "alérgico",
                                "alérgico", "alérgico", "No alérgico",
                                "alérgico", "No alérgico", "No alérgico",
                                "alérgico","alérgico", "alérgico"),
                     mutacion = c(FALSE, FALSE, FALSE, FALSE, TRUE, FALSE,
                                  FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, TRUE,
                                  TRUE, FALSE,  FALSE, TRUE, FALSE, TRUE, FALSE,
                                  TRUE, FALSE,FALSE, FALSE, TRUE, FALSE, FALSE,
                                  TRUE, FALSE, TRUE))

paste("El test de fisher trabaja con tablas")

## [1] "El test de fisher trabaja con tablas"

tabla <- table(datos$sujeto, datos$mutacion, dnn = c("Sujeto", "Estado gen"))
tabla

##              Estado gen
## Sujeto        FALSE TRUE
##   alérgico        6   10
##   No alérgico    11    3

fisher.test(x = tabla, alternative = "two.sided")

##
##  Fisher's Exact Test for Count Data
##
## data:  tabla
## p-value = 0.03293
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.02195148 1.03427479
## sample estimates:
## odds ratio
##  0.1749975

library(vcd)

## Loading required package: grid

paste0("Fuerza de asociación")

## [1] "Fuerza de asociación"

assocstats(x = tabla)

##                     X^2 df P(> X^2)
## Likelihood Ratio 5.3356  1 0.020894
## Pearson          5.1293  1 0.023525
##
## Phi-Coefficient   : 0.413
## Contingency Coeff.: 0.382
## Cramer's V        : 0.413

En este ejemplo no se satisface la condición de frecuencias marginales fijas y por lo tanto el test de Fisher no es exacto. Aun así, hay evidencias para rechazar la H0 y considerar que las dos variables sí están relacionadas. El tamaño de la fuerza de asociación (tamaño de efecto) cuantificado por phi o Cramers V es mediano.

fila1 <- c(81, 103, 147)
fila2 <- c(359, 326, 277)
tabla <- as.table(rbind(fila1, fila2))
dimnames(tabla) = list(Peso = c("Obeso","No obeso"),
                       Estado_civil = c("soltero","pareja","casado"))
tabla

##           Estado_civil
## Peso       soltero pareja casado
##   Obeso         81    103    147
##   No obeso     359    326    277

chisq.test(x = tabla)

##
##  Pearson's Chi-squared test
##
## data:  tabla
## X-squared = 30.829, df = 2, p-value = 2.021e-07

paste("Solución mediante simulación")

## [1] "Solución mediante simulación"

chisq.test(tabla, simulate.p.value = TRUE, B = 5000)

##
##  Pearson's Chi-squared test with simulated p-value (based on 5000
##  replicates)
##
## data:  tabla
## X-squared = 30.829, df = NA, p-value = 2e-04

chisq.test(x = tabla)$residuals

##           Estado_civil
## Peso          soltero     pareja     casado
##   Obeso    -2.9809729 -0.6509186  3.6914422
##   No obeso  1.7485759  0.3818151 -2.1653222

chisq.test(x = tabla)$stdres

##           Estado_civil
## Peso          soltero     pareja     casado
##   Obeso    -4.2549504 -0.9231681  5.2203206
##   No obeso  4.2549504  0.9231681 -5.2203206

solteros_casados <- tabla[, c(1,3)]
solteros_casados

##           Estado_civil
## Peso       soltero casado
##   Obeso         81    147
##   No obeso     359    277

chisq.test(solteros_casados)

##
##  Pearson's Chi-squared test with Yates' continuity correction
##
## data:  solteros_casados
## X-squared = 28.56, df = 1, p-value = 9.083e-08

assocstats(solteros_casados)

##                     X^2 df   P(> X^2)
## Likelihood Ratio 29.687  1 5.0775e-08
## Pearson          29.391  1 5.9140e-08
##
## Phi-Coefficient   : 0.184
## Contingency Coeff.: 0.181
## Cramer's V        : 0.184

Test de McNemar DATOS PAREADOS

El test de McNemar es la alternativa a los test ??2 de Pearson y al test exacto de Fisher cuando: los datos son pareados, se trata de tablas 2x2 y ambas variables son dicotómicas (binomiales). El test de McNemar estudia si la probabilidad de evento verdadero para una variable es igual en los dos niveles de otra variable.

Tabla de McNemar

Supóngase que se quiere comprobar si un tratamiento de hipnosis es capaz de hacer que las personas contesten ???Sí??? con mayor frecuencia. Para ello se selecciona un grupo de individuos a los que se les realiza una pregunta cuya respuesta puede ser SI/NO antes y después de someterse al tratamiento de hipnosis.

Ho: La respuesta de los sujetos es independiente del tratamiento. La proporción de sujetos que pasan de responder Sí a No es igual a la proporción de sujetos que pasan de responder No a Sí.

Ha: El tratamiento sí influye en la respuesta por lo que la proporción de sujetos que pasan de No a Sí es diferente de la de sujetos que pasan de Sí a No.

datos <- data.frame( sujeto = rep(1:15, each = 2),
                     tratamiento = c("pre" , "post" , "pre" , "post" , "pre" ,
                                     "post" , "pre" , "post" , "pre" , "post" ,
                                     "pre" , "post" , "pre" , "post" , "pre" ,
                                     "post" , "pre" , "post" , "pre" , "post" ,
                                     "pre" , "post" , "pre" , "post" , "pre" ,
                                     "post" , "pre" , "post" , "pre" , "post"),
                     respuesta = c("NO", "SI", "SI", "SI", "NO", "SI", "SI",
                                   "NO", "SI", "SI", "NO", "SI", "NO", "SI",
                                   "NO", "SI", "NO", "SI", "SI", "SI", "NO",
                                   "NO", "SI", "SI", "NO", "SI", "NO", "NO",
                                   "NO", "SI"))
library(tidyr)
datos <- spread(data = datos, key = tratamiento, value = respuesta)
head(datos)

##   sujeto post pre
## 1      1   SI  NO
## 2      2   SI  SI
## 3      3   SI  NO
## 4      4   NO  SI
## 5      5   SI  SI
## 6      6   SI  NO

tabla <- table(Pre_Tratamiento = datos$pre, Post_Tratamiento = datos$post)
tabla

##                Post_Tratamiento
## Pre_Tratamiento NO SI
##              NO  2  8
##              SI  1  4

mcnemar.test(tabla)

##
##  McNemar's Chi-squared test with continuity correction
##
## data:  tabla
## McNemar's chi-squared = 4, df = 1, p-value = 0.0455

binom.test(x = 1, n = 1 + 8, p = 0.5 )

##
##  Exact binomial test
##
## data:  1 and 1 + 8
## number of successes = 1, number of trials = 9, p-value = 0.03906
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.002809137 0.482496515
## sample estimates:
## probability of success
##              0.1111111

Existen evidencias significativas para rechazar Ho en favor de que sí existe relación entre el tratamiento y la respuesta de los sujetos.

Respuesta <- c(1,0,1,0,0,1,0,1,0,0,1,1,0,0,1,1,1,1,0,0,1,0,1,1,0,1,1,0,1,1)
Sujeto <- factor(c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7,8,8,8,9,9,9,10,
                   10,10))
Canal <- factor(rep(c("inicial","anuncio","internet"), 10))
datos <- data.frame(Sujeto,Canal,Respuesta)

 require(coin)
symmetry_test(Respuesta ~ factor(Canal) | factor(Sujeto), data = datos,teststat = "quad") 

##
##  Asymptotic General Symmetry Test
##
## data:  Respuesta by
##   factor(Canal) (anuncio, inicial, internet)
##   stratified by factor(Sujeto)
## chi-squared = 8.2222, df = 2, p-value = 0.01639

 require(coin)
symmetry_test(Respuesta ~ factor(Canal) | factor(Sujeto), data = datos,teststat = "quad") 

##
##  Asymptotic General Symmetry Test
##
## data:  Respuesta by
##   factor(Canal) (anuncio, inicial, internet)
##   stratified by factor(Sujeto)
## chi-squared = 8.2222, df = 2, p-value = 0.01639

require(coin)
require(dplyr)

## Loading required package: dplyr

##
## Attaching package: 'dplyr'

## The following object is masked from 'package:car':
##
##     recode

## The following objects are masked from 'package:Hmisc':
##
##     src, summarize

## The following object is masked from 'package:MASS':
##
##     select

## The following objects are masked from 'package:stats':
##
##     filter, lag

## The following objects are masked from 'package:base':
##
##     intersect, setdiff, setequal, union

require(tidyr)
# Inicial vs anuncio
inicial_anuncio <- datos %>% filter(Canal %in% c("inicial","anuncio")) %>% 
                   spread(key = Canal, value = Respuesta)
tabla_inicial_anuncio <- table(inicial = inicial_anuncio$inicial,
                               anuncio = inicial_anuncio$anuncio)
tabla_inicial_anuncio

##        anuncio
## inicial 0 1
##       0 3 5
##       1 1 1

p_value_12 <- mcnemar.test(tabla_inicial_anuncio)
p_value_12

##
##  McNemar's Chi-squared test with continuity correction
##
## data:  tabla_inicial_anuncio
## McNemar's chi-squared = 1.5, df = 1, p-value = 0.2207

# Inicial vs internet
inicial_internet <- datos %>% filter(Canal %in% c("inicial","internet")) %>% 
                    spread( key = Canal, value = Respuesta)
tabla_inicial_internet <- table(inicial = inicial_internet$inicial,
                                internet = inicial_internet$internet)
tabla_inicial_internet

##        internet
## inicial 0 1
##       0 1 7
##       1 0 2

p_value_13 <- mcnemar.test(tabla_inicial_internet)
p_value_13

##
##  McNemar's Chi-squared test with continuity correction
##
## data:  tabla_inicial_internet
## McNemar's chi-squared = 5.1429, df = 1, p-value = 0.02334

# anuncio vs internet
anuncio_internet <- datos %>% filter(Canal %in% c("anuncio","internet")) %>%
                    spread(key = Canal, value = Respuesta)
tabla_anuncio_internet <- table(anuncio = anuncio_internet$anuncio,
                                internet = anuncio_internet$internet)
tabla_anuncio_internet

##        internet
## anuncio 0 1
##       0 0 4
##       1 1 5

p_value_23 <- mcnemar.test(tabla_anuncio_internet)
p_value_23

##
##  McNemar's Chi-squared test with continuity correction
##
## data:  tabla_anuncio_internet
## McNemar's chi-squared = 0.8, df = 1, p-value = 0.3711

p.adjust(c(p_value_12$p.value, p_value_13$p.value, p_value_23$p.value), method = "holm")

## [1] 0.44134272 0.07002661 0.44134272

P-VALUE

Hemos visto cómo calcular comparar las media de dos grupos y decir cuál es su diferencia en un 95% de los casos con un inervalo de confianza, y hemos visto como contrastar los resultados de dos grupos, comporbando si el valor de uno estaba o no dentro del 5% de posibilidades del otro para poder descartar una hipótesis frente a la otra. Estos dos conceptos los unimos para porporcionar en un contraste de hipótesis el p-value, con el que compararmos nuestra hipótesis con la posibilidad de que al azar obtengamos ese mismo valor, proporcionando un intervalo de confianza.

Nos sirve para dictaminar que el resultado obtenido es estadísticamnete significativo.

qownnotes-media-AeDUXI

Es la probabilidad de obtener un resultado al menos tan extremo como el que realmente se ha obtenido (valor del estadístico calculado), suponiendo que la hipótesis nula es cierta. EL pvalue por tanto es una medida de significación estadística (un resultado o efecto es estadísticamente significativo cuando es improbable que haya sido debido al azar).

Una característica de p es que su valor depende del tamaño de la muestra a partir de la cual se determina: a mayor tamaño muestral mayor es la probabilidad de encontrar diferencias significativas entre los grupos que se comparan. De ahí la importancia de que en todo estudio deba realizarse el cálculo del tamaño muestral necesario en base a las caractrísticas de dicho estudio con anterioridad a la recogida y el análisis de los datos.

qownnotes-media-WX4036

Todos los libros de estadística tienen una tabla detrás con los pvalues asociado a cada uno de los valores de Zscore de distintas distribuciones, en este caso distribución normal

qownnotes-media-Tt4036

En primer paso por tanto para probar nuestra hipótesis será crear la H_o para tener una base sobre la que medir una alternativa a nuestra H_a con los datos obervados, calcular un test estadístico desde los datos dados, y finalmente comparar este test estadístico con la distribución hipotética. Esta comparación nos dice cómo de alejado el valor del test está de la hipótesis alternativa. Pvalue es la probabilidad sobre la H_0 de obtener evidencia tanto o más alejada de test estadístico obtenido de los datos observados en la dirección de la hipótesis alternativa.

Por tanto si pvalue es pequeño, o H_o es cierto y los datos observados tienen eventos raros o H_o es falso. Pvalye es el valor más pequeño de alfa en el cual rechazaríamos la H_0

P-value permite elegir un porcentaje de error Typo I , comparando ese porcenaje de error si P-value es menor rechazar la hipotesis nula. Hay muchos mal entendidos en su interpretación, es por eso que tiene muchos defensores y muchos detractores.

Vamos a calcular, cómo de inusal son nuestros dato estimado asumiento que Ho es cierto

qownnotes-media-PC7296

Para calcular pvalue en una distribución nomal, existe la función pnorm() con los distintos parámetros de la distribución, uno de ellos, lower.trail = TRUE/FALSE le decimos si queremos calcularlo por un sólo extremo o por los dos extremos.En el ejemplo queríamos ver el valor de alfa = 2, si sólo cogemos un extremo tendremos alfa >2

qownnotes-media-H11100

Si la ditribución es binomial pbinom.

qownnotes-media-wSCjdp

Si la distribución es de poisson ppois()

qownnotes-media-HOxTYi

qownnotes-media-u11100

Proporcionando el valor de pvalue en lugar de un alfa (probabilidad de cometer error tipo I) permites a los revisores de tu trabajo realizar hipótesis con el alpha que ellos elijan. La norma general será que si pvalue es menor que el alfa, rechazas H_o.

Si la hipótesis abarca ambos extremos, lo que se suele hacer es doblar el pvalue.

qownnotes-media-CPmAdQ

qownnotes-media-ptQNCu

qownnotes-media-dQbvin

qownnotes-media-tgvovm

qownnotes-media-krqBtz

qownnotes-media-QHoCgi

qownnotes-media-sNVxow

qownnotes-media-rMQxye

qownnotes-media-MKVoiM

qownnotes-media-kCHvKO

TEST

qownnotes-media-xOUDhE

qownnotes-media-uOtGOY

qownnotes-media-nLCkjG

qownnotes-media-SaFnex

qownnotes-media-ALBoXw

qownnotes-media-cPFSzZ

qownnotes-media-hbyEwA

PODER DE CONTRASTE

Poder es la probabilidad de rechazar la hipótesis nula cuando es falsa. El error tipo II es llamada Beta, es la proabilidad de no rechazar H_o cuando en realidad es falsa, por tanto el poder = 1-beta.

Se usa por ejemplo, cuando quieres determinar si el tamaño de la muestra es sufientemente grande para arrojar resultados significativos que no se consideren por azar.

En el ejemplo de eventos nocturnos que causan apnea del sueño, la zona azul que queda a la derecha de la línea vertical es el poder de contraste

qownnotes-media-C16176

El poder crecerá cuanto más grande sea la muestra y contra más alejado esté el velor que intentamos inferir del calculado por H_0, en el dibujo de arriba, cuando mayor sea la media de la curva azu, más se aleja del alfa que es la recta y negral mayor área caerá sobre ella a partir de la línea vertical que representa el 5% de H_0 (alfa o error tipo I). En este gráfico vemos cómo según aumenta el valor de media, el poder aumenta, y de manera mucho más incremental si aumentamos el número de muestras

qownnotes-media-sv4036

Cuanto más pequeá sea la desviación, más delgada será la campana, por lo que las curvas se solapan menos y su poder será mayor. osea poder de contraste decrecerá según aumente la desviación estándar ya que significa más variabilidad en los datos y más se solapen H_0 y H_A

qownnotes-media-Y16176

qownnotes-media-zvFTEF

Podemos calcular el poder de dos maneras, o bien manual o bien con la funcion t.test.

qownnotes-media-izEfWo

Para calcularlo de manera manual, usaremos la función pnorm():

z <- qnorm(.95) pnorm(30+z, mean = 32, lower.tail=FALSE)

Esto nos devolverá 0.64, que signiica que con un 64% podemos rechazar H_o cuando es falsa.

Tenemos la función power.t.test que nos realiza esta tarea para distribuciones Tstudent

qownnotes-media-fJ7448

En este caso vemos que tanto Delta (diferencia entre medias) como sd, si se cambian proporcionalmente, el cambio no afecta al poder. Este ration (Delta/sd) se llama effect size.

Se puede utilizar la misma fórmula para saber el número de individuos mínimos necesitamos para cumplir con un power concreto, esto es muy útil para saber el tamaño mínimo de la muestra

qownnotes-media-in7448

EJERCICIOS

Attaching package: ‘swirl’

The following objects are masked from ‘package:shinyjs’:

info, reset

Warning message: package ‘swirl’ was built under R version 3.3.3

swirl()

What shall I call you? miguel

1: Statistical Inference Hypothesis Testing 2: Statistical Inference Probability2 3: Statistical Inference T Confidence Intervals 4: Statistical Inference Variance 5: No. Let me start something new.

Selection: 5

1: Exploratory Data Analysis 2: Getting and Cleaning Data 3: R Programming 4: Statistical Inference 5: Take me to the swirl course repository!

Selection: 4

1: Introduction 2: Probability1 3: Probability2 4: ConditionalProbability 5: Expectations
6: Variance 7: CommonDistros 8: Asymptotics 9: T Confidence Intervals 10: Hypothesis Testing
11: P Values 12: Power 13: Multiple Testing 14: Resampling

Selection: 12

| | 0%

…

|= | 1%

…

|=== | 2%

…

|==== | 3%

…

|====== | 4%

…

|======= | 5%

1: Rejecting a true null hypothesis 2: Miscalculating a t score 3: Misspelling the word hypothesis 4: Accepting a false null hypothesis

Selection: 4

|========= | 7%

…

|========== | 8%

…

|============ | 9%

1: (X’-30)/(s^2/n) 2: 30/(s/sqrt(n)) 3: (X’-30)/(s/sqrt(n)) 4: X’/(s^2/n)

Selection: 3

|============= | 10%

1: a standard measure 2: a standard error 3: a standard sample 4: a standard variance 5: a standard bearer

Selection: 2

|=============== | 11%

…

|================ | 12%

…

|================== | 13%

…

|=================== | 14%

…

|===================== | 15%

…

|====================== | 16%

…

|======================== | 17%

1: the blue 2: the red

Selection: 1

1: the blue 2: the red

Selection: 2

|========================= | 18%

1: the red 2: the blue

Selection: 2

|=========================== | 20%

1: 28 2: 32 3: 36 4: 30

Selection: 2

|============================ | 21%

…

|============================== | 22%

…

|=============================== | 23%

…

|================================= | 24%

…

|================================== | 25%

…

|==================================== | 26%

myplot(34)

|===================================== | 27%

myplot(33.3)

|======================================= | 28%

myplot(30)

|======================================== | 29%

…

|========================================== | 30%

…

|=========================================== | 32%

…

|============================================= | 33%

…

|============================================== | 34%

myplot(28)

|================================================ | 35%

…

|================================================= | 36%

…

|=================================================== | 37%

1: power is independent of mu_a 2: as it gets bigger, it gets less powerful 3: as it gets bigger, it gets more powerful

Selection: 3

|==================================================== | 38%

1: as it gets bigger, it gets more powerful 2: as it gets bigger, it gets less powerful 3: power is independent of sample size

Selection: 1

|====================================================== | 39%

…

|======================================================= | 40%

…

|========================================================= | 41%

…

|========================================================== | 42%

…

|============================================================ | 43%

1: returns the area under the red curve to the left of the line 2: returns the area under the red curve to the right of the line 3: returns the area under the blue curve to the left of the line 4: returns the area under the blue curve to the right of the line

Selection: 4

|============================================================= | 45%

1: an answer dependent on beta 2: an answer less than .50 3: an answer dependent on alpha 4: an answer greater than 1

Selection: 2

|=============================================================== | 46%

z <- qnorm(.95)

|================================================================ | 47%

z <- pnorm(30+z,mean=30,lower.tail = FALSE)

pnorm(30+z,mean=30,lower.tail = FALSE) [1] 0.05

|================================================================== | 48%

pnorm(30+z,mean=32,lower.tail = FALSE) [1] 0.63876

|=================================================================== | 49%

…

|==================================================================== | 50%

…

|====================================================================== | 51%

1: less than 2: the same 3: greater

Selection: 1

|======================================================================= | 52%

pnorm(30+z,mean=32,sd=1,lower.tail = FALSE) [1] 0.63876

|========================================================================= | 53%

pnorm(30+z,mean=32,sd=2,lower.tail = FALSE) [1] 0.5704709

pnorm(30+z*2,mean=32,sd=2,lower.tail = FALSE) [1] 0.259511

|========================================================================== | 54%

1: it doesn’t change 2: it decreases 3: it increases

Selection: 3

|============================================================================ | 55%

…

|============================================================================= | 57%

1: it decreases 2: it increases 3: it doesn’t change

Selection: 2

|=============================================================================== | 58%

…

|================================================================================ | 59%

1: it doesn’t change 2: it increases 3: it decreases

Selection: 3

|================================================================================== | 60%

…

|=================================================================================== | 61%

1: it decreases 2: it doesn’t change 3: it increases

Selection: 3

|===================================================================================== | 62%

…

|====================================================================================== | 63%

1: less than 2: greater 3: they’re the same

Selection: 3

1: less than 2: greater 3: they’re the same

Selection: 2

|======================================================================================== | 64%

1: No 2: Yes

Selection: 2

|========================================================================================= | 65%

…

|=========================================================================================== | 66%

1: alpha 2: mu_0 3: beta 4: mu_a

Selection: 2

|============================================================================================ | 67%

1: mu_0 2: alpha 3: mu_a 4: beta

Selection: 2

|============================================================================================== | 68%

…

|=============================================================================================== | 70%

…

|================================================================================================= | 71%

…

|================================================================================================== | 72%

…

|==================================================================================================== | 73%

…

|===================================================================================================== | 74%

…

|======================================================================================================= | 75%

…

|======================================================================================================== | 76%

power.t.test(n = 16, delta = 2 / 4, sd=1, type = “one.sample”,alt = “one.sided”)$power [1] 0.6040329

|========================================================================================================== | 77%

power.t.test(n = 16, delta = 2, sd=4, type = “one.sample”,alt = “one.sided”)$power [1] 0.6040329

|=========================================================================================================== | 78%

power.t.test(n = 16, delta = 100, sd=200, type = “one.sample”,alt = “one.sided”)$power [1] 0.6040329

|============================================================================================================= | 79%

…

|============================================================================================================== | 80%

power.t.test(power = .8, delta = 2 / 4, sd=1, type = “one.sample”, alt = “one.sided”)$n [1] 26.13751

|================================================================================================================ | 82%

power.t.test(power = .8, delta = 2, sd=4, type = “one.sample”, alt = “one.sided”)$n [1] 26.13751

|================================================================================================================= | 83%

power.t.test(power = .8, delta = 100, sd=200, type = “one.sample”, alt = “one.sided”)$n [1] 26.13751

|=================================================================================================================== | 84%

power.t.test(power = .8, n = 26, sd=1, type = “one.sample”, alt = “one.sided”)

 One-sample t test power calculation

          n = 26
      delta = 0.5013986
         sd = 1
  sig.level = 0.05
      power = 0.8
alternative = one.sided

power.t.test(power = .8, n = 26, sd=1, type = “one.sample”, alt = “one.sided”)$delta [1] 0.5013986

|==================================================================================================================== | 85%

power.t.test(power = .8, n = 27, sd=1, type = “one.sample”, alt = “one.sided”)$delta [1] 0.4914855

|====================================================================================================================== | 86%

1: delta will halve 2: delta will double 3: delta won’t change

Selection: 1

1: delta won’t change 2: delta will halve 3: delta will double

Selection: 1

1: delta will halve 2: delta won’t change 3: delta will double

Selection: 3

|======================================================================================================================= | 87%

…

|========================================================================================================================= | 88%

1: gamma 2: None of the others 3: alpha 4: beta 5: delta

Selection: 5

1: gamma 2: delta 3: None of the others 4: alpha 5: beta

Selection: 4

|========================================================================================================================== | 89%

1: accepting a true hypothesis 2: rejecting a false hypothesis 3: rejecting a true hypothesis 4: accepting a false hypothesis

Selection: 4

|============================================================================================================================ | 90%

1: delta 2: thrilling 3: None of the others 4: alpha 5: beta 6: gamma

Selection: 3

|============================================================================================================================= | 91%

1: No 2: Yes

Selection: 1

|=============================================================================================================================== | 92%

1: huh? 2: mu_a-mu_0 < 0 3: mu_a-mu_0=0 4: mu_0-mu_a > 0

Selection: 3

|================================================================================================================================ | 93%

1: mu_0=mu_a maximizes the power 2: the smaller mu_0-mu_a the more powerful the test 3: the smaller mu_a-mu_0 the more powerful the test

Selection: 2

1: the smaller mu_0-mu_a the more powerful the test 2: mu_0=mu_a maximizes the power 3: the smaller mu_a-mu_0 the more powerful the test

Selection: 3

|================================================================================================================================== | 95%

1: (mu_a - mu_0) / sqrt(sigma) 2: (mu_a - mu_0) / sqrt(n) 3: (mu_a - mu_0) / n 4: (mu_a - mu_0) / sigma

Selection: 4

|=================================================================================================================================== | 96%

1: False 2: True

Selection: 2

|===================================================================================================================================== | 97%

1: True 2: False

Selection: 2

|====================================================================================================================================== | 98%

1: False 2: True

Selection: 1

|======================================================================================================================================== | 99%

…

|=========================================================================================================================================| 100%

1: No 2: Yes

COMPARACIONES MÚLTIPLES (MULTIPLE TESTING)

Es un todo un subcampo detro de la inferencia estadística. Trata de reducir y ajustar los errores tipo I y II para evitar falsos positivos (rechazar H_o cuando no se debe) y falsos negativos (rechazar H_a cuando no se debe).

Imaginemos que queremos probar una hipótesis sobre múltiples grupos. Lo que haremos será comparar las medias obtenidas de todos esos grupos.

Un ejemplo claro que podemos tener en esta situación es que, pongamos que hacemos 10.000 test, si calculamos un PValue de manera que P< alfa estaremos controlando los falsos positivos en un nivel de media alfa. Suponemos P< 0.05 significante, por lo que 10000*0.05 = 500 posibles falsos positivos, que parece un número muy alto. Cómo evitamos esto?

Vamos a jugar con ALFA equilibrando la reducción del número de falsos positivos con que tengamos una ALFA muy baja y no podamos rechazar nunca Ho

qownnotes-media-y17684

FWER Para controlar el falso positivo: Family wise error rate (FWER) que es la probabilidad de encontrar al menos 1 falso positivo. Se puede intentar controlar este FWER con la corrección de Bonferroni. Lo bueno de este método es que es muy rápido, pero baja tanto ALFA que nos será muy dificil rechazar Ho

qownnotes-media-X17684

FDR

Otro método para controlar los falsos positivos es False discovery rate (FDR) most popular in genomics E(V/R) = Ratio de falsos positivos (número de falsos positivos divido número de positivos). En este caso se aplica la corrección BH (Benjamini-Hochberg). Es fácil de calcular y aunque permite más falsos positivos es menos conservador y permite rechazar más facilemnte Ho

qownnotes-media-Z17684

En esta gráfica, se ve que el número de resultados significantes sin corrección son 4, con la corrección HC baja a 3 y con la de bonferoni a 2

qownnotes-media-F17684

Adjusted Pvalues

Hay otra tercera aproximación que es en lugar de buscar un umbral a ALFA ajustar pvalues de las muestras

qownnotes-media-Y17684

qownnotes-media-fMlJWx

qownnotes-media-KNlCiF

qownnotes-media-WpNVqa

qownnotes-media-cGLCMd

qownnotes-media-ubeOBL

qownnotes-media-RrJxGj

En R hay una función llamada p.adjust que hace esto, teniendo un vector de mil pvalues:

EJERCICIOS

…

|===== | 5%

…

|======= | 7%

…

|========= | 8%

…

|========== | 10%

1: rejecting a false hypothesis 2: failing to reject a true hypothesis 3: rejecting a true hypothesis 4: failing to reject a false hypothesis

Selection: 4

1: rejecting a true hypothesis 2: failing to reject a true hypothesis 3: rejecting a false hypothesis 4: failing to reject a false hypothesis

Selection: 3

1: rejecting a true hypothesis 2: failing to reject a false hypothesis 3: failing to reject a true hypothesis 4: rejecting a false hypothesis

Selection: 1

|============ | 11%

1: convicting an innocent person 2: letting the indicted off on a technicality 3: acquitting a guilty person

Selection: 3

1: letting the indicted off on a technicality 2: acquitting a guilty person 3: convicting an innocent person

Selection: 3

|============== | 13%

1: rejecting a false hypothesis 2: failing to reject a false hypothesis 3: failing to reject a true hypothesis 4: rejecting a true hypothesis

Selection: 1

1: failing to reject a false hypothesis 2: rejecting a true hypothesis 3: rejecting a false hypothesis 4: failing to reject a true hypothesis

Selection: 1

|=============== | 15%

1: letting the indicted off on a technicality 2: convicting an innocent person 3: acquitting a guilty person

Selection: 3

|================= | 16%

1: tells us the origins of the number 0 2: represents the status_quo and is assumed true 3: is a big nothing that statisticians like to gossip about 4: is never true

Selection: 2

|=================== | 18%

…

|==================== | 20%

1: revised 2: renamed the aleph null hypothesis 3: rejected 4: accepted

Selection: 3

|====================== | 21%

…

|======================== | 23%

1: A,B,C 2: S,T,U,V 3: m_0, and m 4: m and R

Selection: 4

|========================== | 25%

1: Type I 2: Type III 3: Type II 4: a serious one

Selection: 1

|=========================== | 26%

…

|============================= | 28%

1: Type III 2: Type II 3: a serious one 4: Type I

Selection: 2

|=============================== | 30%

…

|================================ | 31%

…

|================================== | 33%

…

|==================================== | 34%

…

|====================================== | 36%

1: false alarm rate 2: the Type II rate 3: a thorn 4: a rose

Selection: 1

|======================================= | 38%

1: the Type I error rate 2: the Type II error rate 3: a thorny rose

Selection: 1

|========================================= | 39%

…

|=========================================== | 41%

…

|============================================ | 43%

…

|============================================== | 44%

1: 50000 2: 5000 3: 50 4: 500

Selection: 4

|================================================ | 46%

…

|================================================= | 48%

…

|=================================================== | 49%

…

|===================================================== | 51%

1: too many results will fail 2: too many results will pass 3: requires too much math

Selection: 1

|======================================================= | 52%

…

|======================================================== | 54%

…

|========================================================== | 56%

…

|============================================================ | 57%

…

|============================================================= | 59%

…

|=============================================================== | 61%

1: 8 2: 6 3: 2 4: 4

Selection: 4

|================================================================= | 62%

1: 6 2: 2 3: 4 4: 8

Selection: 2

|================================================================== | 64%

1: 1 2: 7 3: 5 4: 3

Selection: 4

|==================================================================== | 66%

…

|====================================================================== | 67%

…

|======================================================================== | 69%

…

|========================================================================= | 70%

…

|=========================================================================== | 72%

head(pValues) [1] 0.5334915 0.2765785 0.8380943 0.6721730 0.8122037 0.4078675

|============================================================================= | 74%

sum() [1] 0

sum(pValues < 0.05) [1] 51

|============================================================================== | 75%

p.adjust(pValue,method=“bonferroni”) Error in p.adjust(pValue, method = “bonferroni”) : object ‘pValue’ not found p.adjust(pValues,method=“bonferroni”)

sum(p.adjust(pValues,method=“bonferroni”) < 0.05) [1] 0

|================================================================================ | 77%

sum(p.adjust(pValues,method=“BH”) < 0.05) [1] 0

|================================================================================== | 79%

…

|==================================================================================== | 80%

tail(trueStatus) [1] “not zero” “not zero” “not zero” “not zero” “not zero” “not zero”

|===================================================================================== | 82%

table(trueStatus) trueStatus not zero zero 500 500

table(pValues2 < 0.05, trueStatus) trueStatus not zero zero FALSE 0 476 TRUE 500 24

|======================================================================================= | 84%

…

|========================================================================================= | 85%

4.8 [1] 4.8

24/500 [1] 0.048

|========================================================================================== | 87%

…

|============================================================================================ | 89%

table(pValues2 < 0.05, trueStatus) trueStatus not zero zero FALSE 0 476 TRUE 500 24

table(p.adjust(pValues2,method=“bonferroni”) < 0.05, trueStatus) trueStatus not zero zero FALSE 23 500 TRUE 477 0

|============================================================================================== | 90%

…

|=============================================================================================== | 92%

table(p.adjust(pValues2,method=“BH”) < 0.05, trueStatus) trueStatus not zero zero FALSE 0 487 TRUE 500 13

|================================================================================================= | 93%

…

|=================================================================================================== | 95%

…

|===================================================================================================== | 97%

…

|====================================================================================================== | 98%

…

|========================================================================================================| 100%

RESAMPLING INFERENCE

Hasta ahora hemos visto métodos de inferencia estadística para obtener los parámetros de una población en base a su distrubución, pero ¿qué pasa si desconocemos su distribución?¿Cómo podemos obterner la media y varianza poblacional por ejemplo al tirar un dado si sabemos que está trucado? De lo que se trata entonces es de realizar inferencias sobre nuestros datos tomando muestras de nuestros propios datos, es decir haciendo múltiples selecciones aleatorias sobre los datos que hemos recogido, para poder obterner la media de sus medias.

BOOTSTRAPING

Es una técnica que permite simplificar la obtención de los intervalos de confianza en lugar de con métodos matemáticos, calcularlos a través de múltiples muestras aleatorias tomadas de las muestras sin sacar del espacio muestral en cada iteración la muesta obtenida (reemplazo).

Es una técnica no paramétrica y nos serviría para poder detallar en nuestro estudio un intervalo de confianza y el error estándar sobre una única medición tomada, iterando una y otra vez sobre esa muesta para obtener de manera aleatoria uno se sus valores y así simular la población.

qownnotes-media-JkhSIH

qownnotes-media-tanpiX

qownnotes-media-oTTWHa

qownnotes-media-HmdAlZ

qownnotes-media-MvzZKg

qownnotes-media-ugudYv

qownnotes-media-CGFVKT

qownnotes-media-zblFDp

PERMUTATION TEST

Esta es otra técnica de inferencia por muestreo, en la que en lugar de utilizar T-Student o Chicuadrado, como no sabemos la distribución de los datos comparamos los grupos haciendo permutaciones aleatorias entre ellos mismos y lo comparamos con los resultados más extremos que podemos obtener de su comparativa inicial.

A grandes rasgos, se restan las medias de ambos grupos y se comprueba con otro grupo de permutaciones realizadas cuántas restas de las medias de las permutaciones aleatorias supera la medida inicial, y así saber cuán extremo es la diferencia.

qownnotes-media-WVBDaB

qownnotes-media-lrLfbM

qownnotes-media-RRTRbe

qownnotes-media-BxnNue

qownnotes-media-txJbfA

EJERCICIOS

…

|==== | 4%

…

|====== | 6%

…

|======= | 7%

1: a statistic 2: a hypothesis 3: a population 4: observations

Selection: 1

1: a statistic 2: a population 3: observations 4: a hypothesis

Selection: 2

|========= | 8%

1: a statistic 2: a hypothesis 3: observations 4: a population

Selection: 3

|========== | 10%

…

|============ | 11%

1: an entirely new sample 2: a worse sample 3: the original sample permuted 4: a better sample

Selection: 3

|============= | 12%

…

|============== | 14%

3 [1] 3

sum(1:6)/6 [1] 3.5

|================ | 15%

…

|================= | 17%

…

|=================== | 18%

…

|==================== | 19%

…

|====================== | 21%

print(g2)

|======================= | 22%

…

|========================= | 24%

…

|========================== | 25%

1: the most frequent outcome 2: 50th percentile 3: a person who talks to spirits 4: a point halfway between rare and well-done

Selection: 1

1: the most frequent outcome 2: a person who talks to spirits 3: 50th percentile 4: a point halfway between rare and well-done

Selection: 3

|=========================== | 26%

head(sh) [1] 59.77827 63.21404 63.34242 62.79238 64.28113 64.24221

|============================= | 28%

length(sh) [1] 1078

nh [1] 1078

|============================== | 29%

…

|================================ | 31%

…

|================================= | 32%

…

|=================================== | 33%

mean(nh) [1] 1078

median(resampledMedians) [1] 68.61273

|==================================== | 35%

warning messages from top-level task callback ‘mini’ Warning message: In rm(list = c(vName), envir = globalenv()) : object ‘resamples’ not found > mean(sh) [1] 68.68407

median(sh) [1] 68.61582

|====================================== | 36%

…

|======================================= | 38%

…

|======================================== | 39%

…

|========================================== | 40%

…

|=========================================== | 42%

…

|============================================= | 43%

sam <-sample(fh,nh*B,replace = TRUE)

|============================================== | 44%

resam <- matrix

resam <- matrix(sam,B,nh)

|================================================ | 46%

meds <- apply(resam,1,median)

|================================================= | 47%

| Now look at the difference between the median of fh and the median of meds. > median(meds)-median(fh + ) [1] -0.004105

|=================================================== | 49%

sd(meds) [1] 0.1041539

|==================================================== | 50%

sd(resampledMedians) [1] 0.08293445

|===================================================== | 51%

quantile(resampledMedians,c(.025,.975)) 2.5% 97.5% 68.43733 68.80718

|======================================================= | 53%

quantile(meds,c(.025,.975)) 2.5% 97.5% 67.54999 67.94110

|======================================================== | 54%

…

|========================================================== | 56%

…

|=========================================================== | 57%

…

|============================================================= | 58%

…

|============================================================== | 60%

…

|================================================================ | 61%

…

|================================================================= | 62%

…

|================================================================== | 64%

dim(InsectSprays) [1] 72 2

|==================================================================== | 65%

names(InsectSprays) [1] “count” “spray”

|===================================================================== | 67%

…

|======================================================================= | 68%

range(Bdata$count) [1] 7 21

|======================================================================== | 69%

range(Cdata$count) [1] 0 7

|========================================================================== | 71%

…

|=========================================================================== | 72%

BCcounts [1] 11 17 21 11 16 14 17 17 19 21 7 13 0 1 7 2 3 1 2 1 3 0 1 4

|============================================================================= | 74%

group [1] “B” “B” “B” “B” “B” “B” “B” “B” “B” “B” “B” “B” “C” “C” “C” “C” “C” “C” “C” “C” “C” “C” “C” “C”

|============================================================================== | 75%

testStat() Error in mean(w[g == “B”]) : argument “w” is missing, with no default testStat function(w, g) mean(w[g == “B”]) - mean(w[g == “C”]) <environment: 0x0000000021419b48>

|=============================================================================== | 76%

obs <- testStat((BCcounts,group)) Error: unexpected ‘,’ in “obs <- testStat((BCcounts,” obs <- testStat(BCcounts,group)

|================================================================================= | 78%

obs [1] 13.25

|================================================================================== | 79%

apply(Bdata\(count-Cdata\)count,1,mean) Error in apply(Bdata\(count - Cdata\)count, 1, mean) : dim(X) must have a positive length apply(Bdata\(count-Cdata\)count,mean) Error in match.fun(FUN) : argument “FUN” is missing, with no default mean(Bdata\(count-Cdata\)count) [1] 13.25

|==================================================================================== | 81%

1: mean is linear 2: the data is special 3: mathemagic

Selection: 1

|===================================================================================== | 82%

…

|======================================================================================= | 83%

…

|======================================================================================== | 85%

sample(gruop) Error in sample(gruop) : object ‘gruop’ not found sample(group) [1] “C” “C” “B” “C” “C” “B” “C” “B” “C” “B” “B” “B” “C” “B” “B” “C” “C” “B” “C” “C” “B” “C” “B” “B”

|========================================================================================== | 86%

…

|=========================================================================================== | 88%

perms <- sapply(1 : 10000, function(i) testStat(BCcounts,sample(group)))

|============================================================================================ | 89%

mean(pems > pems) Error in mean(pems > pems) : object ‘pems’ not found mean(perms > obs) [1] 0

|============================================================================================== | 90%

…

|=============================================================================================== | 92%

…

|================================================================================================= | 93%

…

|================================================================================================== | 94%

testStat(DEcounts,group) [1] 1.416667

|==================================================================================================== | 96%

perms <- sapply(1 : 10000, function(i) testStat(DEcounts, sample(group)))

|===================================================================================================== | 97%

…

|======================================================================================================= | 99%

TEST

qownnotes-media-MNbivq

qownnotes-media-fFZicc

qownnotes-media-VMctaN

qownnotes-media-DwoZEa

qownnotes-media-oXujCs

qownnotes-media-pPrnhV

qownnotes-media-GcRabm

qownnotes-media-ETEsxE

qownnotes-media-RnXkAC